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' Abstract. We present ROBO, a model and its companion code for the study of the interstellar medium (ISM). 

r- { , The aim is to provide an accurate description of the physical evolution of the ISM and to set the ground for 

an ancillary tool to be inserted in NBody-Tree-SPH (NB-TSPH) simulations of large-scale structures in the 
cosmological context or of the formation and evolution of individual galaxies. The ISM model consists of gas and 
dust. The gas chemical composition is regulated by a network of reactions that includes a large number of species 
(hydrogen and deuterium-based molecules, helium, and metals). New reaction rates for the charge transfer in 
and H2 collisions are presented. The dust contains the standard mixture of carbonaceous grains (graphite grains 
and PAHs) and silicates. In our model dust are formed and destroyed by several processes. The model accurately 
treats the cooling process, based on several physical mechanisms, and cooling functions recently reported in the 
literature. The model is applied to a wide range of the input parameters and the results for important quantities 
i-C ■ describing the physical state of the gas and dust are presented. The results are organized in a database suited 

to the artificial neural networks (ANNs). Once trained, the ANNs yield the same results obtained by ROBO 
Q , with great accuracy. We plan to develop ANNs suitably tailored for applications to NB-TSPH simulations of 

5— ( ' cosmological structures and/or galaxies. 

■ Key words. ISM: evolution - ISM: dust - galaxies: formation and evolution - methods: numerical 

J> ' 1. Introduction reduce the computational performances of any numerical 
, algorithm (code) that one may adopt to this purpose. This 
Modeling the gas chemistry is an important step towards requires a strategy for optimizing the chemical accuracy 
correctly describing the growth of cosmological structures, of the ISM model and the computational speed, 
the formation and evolution of galaxies, and star forma- In this paper we present a new model of the ISM and 
^ tion in general. For instance, the molecular hydrogen is one the associated code we have developed to explore the ISM 
^ of the most efficient coolants, and its abundance eventu- properties over wide ranges of the physical parameters 
y—t , ally determines the total amount of stars in the Universe, and, at the same time, to cope with the above difficul- 
' Structure growth and galaxy formation and evolution are ties. The model and companion code are named ROBG0. 
customarily investigated by means of large numerical sim- The model deals with an ideal ISM element of unit 
ulations in which a wide set of chemical reactions taking volume, containing gas and dust in arbitrary initial pro- 
place in the ISM should be considered to get and follow portions, whose initial physical conditions are specified 
the key molecules (elemental species in general) eventually by a set of parameters, which is allowed to evolve for a 
governing the efficiency of the star formation and gas cool- given time interval. The history leading the element to 
ing. However, we must face the growing standard complex- that particular initial physical state is not of interest here, 
ity of a typical NB-TSPH model that includes particles The ISM element is mechanically isolated from the host 
of dark matter, particles of baryonic matter (this in the environment; i.e. it does not expand or contract under 
form of stars and gas, in turn divided into several thermal the action of large-scale forces. It can, however, be inter- 
and chemical phases, such as (i) cold, warm and hot, (ii) ested by the passage of shock waves caused by physical 
atomic and molecular, (iii) neutral and ionized), sources phenomena taking place elsewhere (e.g. supernova explo- 
of energy heating and cooling, energy feedback, and eas- sions). Furthermore, it neither acquires nor loses material, 
ily many other physical processes. For this reason, too, a 

detailed chemical description of the ISM would drastically ^ The name means "thing" in some northern Italian dialects. 
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SO the conservation of total mass applies, even if its chem- 
ical composition can change with time. It is immersed in a 
bath of UV radiation generated either by nearby or inter- 
nal stellar sources and in a field of cosmic ray radiation. It 
can generate its own radiation field by internal processes 
and so it has its own temperature, density, and pressure, 
each related by an Equation of State (EoS). If observed 
from outside, it would radiate with a certain spectral en- 
ergy distribution. For the aims of this study, we do not 
need to know the whole spectral energy distribution of 
the radiation field pervading the element, but only its UV 
component. Given these hypotheses and the initial con- 
ditions, the ISM element evolves toward another physical 
state under the action of the internal network of chemical 
reactions changing the relative abundances of the elemen- 
tal species and molecules, the internal heating and cooling 
processes, the UV radiation field, the field of cosmic rays 
and the passage of shock waves. In view of the future appli- 
cations of this model in dynamical simulations of galaxies, 
the integration time interval is chosen in such a way that 
(i) it is long enough to secure that the secular evolution of 
the gas properties is achieved, and (ii) it is short enough to 
secure that the physical properties of the ISM at each in- 
stant are nearly independent of the external variations in 
the large-scale properties of the system hosting the ISM 
element. In other words, ROBO associates a final state 
(another point in the same space) to any initial state (a 
point in the multidimensional space of the physical pa- 
rameters) through a path. The model is like an operator 
determining the vector field of the local transformations 
of the ISM from an initial state to a final one in the space 
of the physical parameters. This is the greatest merit of 
this approach, which secures the wide applicability of the 
model. 

The new ISM model stands on recipes of the inter- 
nal phY^icaljDrocesses fallin g in between thos e developed 
bv lGlover fc Jappsenl (|2007t ) and lSmith et~all (fioOS) . The 
first one follows the thermal and chemical evolution of 
the low-metallicity gas in large numerical simulations. The 
chemical network includes a detailed treatment of H and 
He but neglects the molecules formed with heavy elements 
as the CO molecule. Dust is included to compute its con- 
tribution to the formation of molecular hydrogen, but its 
evolution is not calculated. The chemical code is an on- 
the-fly routine, running as part of a wider code following 
cosmological simulations of structure growth. 



The model proposed by ISmith et all (l2008h uses the 



non-equilibrium treatment for hydrogen-like species and 
the standard equilibrium approximation for all the re- 
maining chemical species. It does not take any type of dust 



into a ccount. To calcul ate the cooling rates , ISmith et al. 



( 2008h use CLOUDY (jPerland et al.l llQQSh and get the 
cooling rates as a function of temperature, density, and 
metallicity. By doing this, it is possible to include a large 
chemical network and a wide set of coolants, but the price 
to pay is that several oversimplifications of the problem are 
mandatory, e.g. the assumption of ioniz ation equilibrium. 
A similar approach has been proposed bv lHocuk fc Soaana 



(|2010t ) using the iMeiierink fc SpaansI (|2005[ ) PDR model 
instead of CLOUDY. 

Our ISM model and associated code ROBO can not 
only describe the gas evolution in great detail but also in- 
cludes large chemical networks and the presence of various 
types of dust which follow the chemistry and the com- 
plex interplay between grain destruction and formation. 
Many gas and dust components are taken into account, 
among which w e recall the molecular hydrogen and the 
metal coolants (ISantoro fc Shi^ l2006t iMaio et all l2007l) 
or HD (jMcGreer fc BrvanlbOoiT To track the evolution of 
these components, the ISM model and ROBO take var- 
ious physical processes into account that may affect the 
behavior of the whole system. For i nstance, dust is ver y 
efficient in forming both H2 and HD (jCazaux et al.ll2008r ). 
Including dust formation and destruction is a formidable 
task. Among other processes, dust can be destroyed by 
shocks that deserve special care to be properly modeled. 
Here we have considered two different approaches. The 
first one makes use of a mean shock speed that is assumed 
to be the same for all the gas particles, neglecting the 
effects of the environment. In reality this is too crude a 
description. The second approach starts from the notions 
that the shocks develop when the motion of the gas parti- 
cles becomes turbulent and that the distribution of turbu- 
lent velocities obeys the one predicted by the Kolmogorov 
law. This is significantly better than the previous case, so 
it is the approach we prefer. 

Finally, we would like to include the gas model (and 
results of ROBO) in numerical simulations of galaxies in 
a simple way. We suppose that a numerical code calcu- 
lates the formation and evolution of a model galaxy from 
an initial stage at high redshift using the standard NB- 
TSPH technique. At each time step, it requires an update 
of the chemical status of the gas particles. This is done for 
every gas particle of the simulation (typically from 10'* up 
to 10'' particles, depending on the simulation under con- 
sideration). We have two methods to our disposal. The 
first one is a real-time chemical updater. This approach 
must use simple physics and a powerful computer in order 
to save computing time. The second one is to use model 
grids, calculated in advance for a wide range of the in- 
put parameters, in such a way as to cover the plausible 
space of the initial conditions. Since increasing the pa- 
rameters also increases the space dimensions of the grids, 
thus making data handling cumbersome, we make use of 
the ANNs technique to get rid of this difficulty. Once the 
ANNs are instructed to reproduce the (ROBO) results as a 
function of the parameters, they should replace ROBO in 
the NB-TSPH simulator o f galaxies. In our case the NB- 
TS PH model is EvoL by iMerlin fc Chlosil (|2nnd. 1200/1) 
and Merlin et al. The use of ANNs can greatly im- 

prove upon this point of difficulty. Here we briefly touch 
upon this problem and l eave the detailed d iscussion of it 
to a forthcoming paper ( Grassi et al.l 2011 ) . 

The plan of the paper is as follows. In Section 2 we 
give a detailed description of the physics behind the ISM 
model and ROBO. This section is divided in three parts: 
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Table 1. Correspondence between elemental species or 
free particles and the indices of the differential equations 
governing the reaction network. 
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the chemistry of the gas phase, the presence of dust grains 
of different types (including their formation by accretion 
and destruction) , and finally the heating and cooling pro- 
cesses. In Section 3 we describe the general characteristics 
of the code. Section H] is dedicated to describing the re- 
sults of the calculations run to validate ROBO. Section 
[5] describes how to include the results of ROBO in the 
NB-TSPH simulations. Finally, some concluding remarks 
are presented in Section [S) 

2. Physical model of the ISM 

In this section we describe the physics of the ISM model 
that is used by ROBO. The problem and associated code 
are divided into three parts, mirrored by the structure of 
this section: gas chemistry, dust, and cooling. All these 
aspects are mutually coupled and overlapped. 

2.1. The chemical network for the ISM 

The kernel of the numerical code modelling the gas chem- 
istry is the network of chemical reactions among different 
elemental species (atoms and molecules, neutral or ion- 
ized dust grains of different types) and free particles (such 
as the electrons), the network of photochemical reactions 
between the above elemental species, and the radiation 
field. In addition to this, we must include all the processes 
creating and destroying molecules and dust grains with 
particular attention to those that are most efficient for 
cooling the gas. For this reason our chemical network fol- 
lows species like H2, HD, and metals (C, O, Si, Fe, and 
their ions). 

The number of reactions depends on how many species 
are tracked and how many details are included in their 
descripti on. A nuinber of codes study the behavior of 
the ISM jCe n '1992; 'Katz et al.lll9 96': 'Calh fc Pallalll998l: 
Anninos et al. 1997: Glover fc Savin 2009 '). each of which 
has a different number of species to follow and a different 
degree of sophistication for the physical processes taken 
into consideration. 



We keep track of the following 27 elemental species or 
molecules plus the free electrons: H, 11+ , H^, H2, H^, D, 
D+, D-, D2, HD, HD+, He, He+, He++, C, C+, CH, CH2, 
CH^, CH;^, CO, O, 0+, Si, Si+, Fe, Fe+, and e". 

These species are divided into four groups. The first 
one contains hydrogen-only based species. The second 
group is composed of deuterium-based species, in which 
HD plays the key role in the gas cooling. The third group 
lists helium and its ions. Carbon, oxygen, silicon, and iron 
with their ions and compounds form the fourth group. The 
free electrons link all the four groups together. The species 
from CH through CH^ are introduced to follow the for- 
mation and destruction of CO to be described below. 

The reactions in which all the above species are in- 
volved are 

— collisional ionization (A + c~ A+ + 2c~), 

— photo-recombination (A+ + ^ A + 7), 

— dissociative recombination (Aj + — > 2A), 

— charge transfer (A+ -f B ^ A + B+), 

— radiative attachment (A + ^ A 

— dissociative attachment (A 

— collisional detachment (A^ 

— mutual neutralization (A+ 

— isotopic exchange (Aj -I- B 

— dissociations by cosmic rays (AB + CR ^ A + B), 

— neutral- neutral (AB -I- AB -)■ A2 -I- B2), 

— ion- neutral (AB"*" -I- AB ABJ -|- A), 

— colhder (AB + C ^ A + B + C), 

— ionizations by field photons (A + 7 — > A+ + c^) 



B- ^ 
c^ 

B- ^ 
■ AB+ 



+ 7), 

AB + e-), 
A + 2e-), 
A + B), 
+ A) 



where A and B are two generic atoms. To these reactions 
we must add those regulating the abundance of CO to be 
described below. 

The chemical network governing the ISM model is a 
classical system of differential equations in which each 
equation is a Cauchy problem of the form 



dni{t) 
dt 



= J2 Rim{T)ni{t)nm(t) - R^AT)n^{t)n■i{t) , (1) 



where ni{t) is the number density of the i-th species 
with known initial value 71^(0). In eqn. ([T|) Rim{T) is 
the rate of the reaction between the l-th and the m-th 
species expressed in units of cm'^/s. Eqn. ([1} is written 
as the sum of all the reactions forming the i-th species 
(^im Rini{T)ni{t)n„i{t)) minus the sum of all the reac- 
tions destroying the i-th species i^jRij{T)ni{t)n,j{t)). 
The indices and m run from 1 to 28 as in the list 

of elemental species and free electrons already mentioned. 
The one-to-one correspondence between the indices in the 
eqn. ([T]) and elemental species is given in Table [T] 

To clarify the meaning of eqn. (jlj we show the case of 
a two-reaction system 



Ao -I- Al A2 + A3 
A4 A5 ^ Ao A3 , 



(2) 
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Table 2. The reacti on rates among hydrogen, deuterium and helium. References: (a) Glover fc Savin ( 2009f ): (b) 



Glover fc Abell (|2008l ). More details within the text. 
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where Ai is the number density of the generic species in 
units of cm"'^. Looking at the density variation of the 
species Aq over the time step dt the density tiaq changes 
as 

duAoit + dt) ^ - koiiT)nA„{t)nAi{t)dt 

+ ki5{T)nA,{t)nA,{t)<it . (3) 

In this equation the positive term is due to the second re- 
action (between A4 and A5) that increases the total quan- 
tity of Aq (and also A3), while the negative term is due 
to the first reaction in which Aq is destroyed because it 
reacts with Ai to form A2 and A3. 

We adopt here a minimal description containing only 
the 28 species and/or the molecular compounds listed 



in Table [U the 64 reactions listed in Tables [5] and [3] 
(where reactions among hydrogen, deuterium, and he- 
lium are considered) , in Table S] (where metals and cos- 
mic rays are involved), in Table [5] (where the reactions 
of the CO mini-network are listed), and finally, the 12 
photochemical processes listed in Table [SI In all the ta- 
bles the gas temperature T is in Kelvin; the electron tem- 
perature Te is in eV. In Table S] the cosmic ray field is 
CcR and it is the rate of H2 ionization by cosmic 
rays. For Table [5] we have T300 = T/SOO. The expres- 
sions in Table [5] are ^ = E/Eth, ^ph{af,yw,Xf,ya,P) = 

[{Xf - 1)2 + 2/^]y0.5P-5.5[i + ^/^]-P, y = ^x)+a}, 

Xf{a, h) = E/a - 6, and xs^ = 1.672 x 10"^. 
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Table 3. Continuation of Tabled References: (a) IClover fc Siwi3 (|2009[) : (b) IClover fc Abel (|2008f ): (c) lAbel et al 
( 19971) . 
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The various types of reaction are described below in 
some detail. Tables [U [31 HI and [5] also sample the chemi- 
cal reactions according to the same four groups in which 
we separated the elemental species depending on the type 
of reaction and the information at our disposal to de- 
rive the reactions rates. In doing this, the notation in use 
may appear rather complex. This is done on purpose, and 
the main motivation is that we want to keep the same 
formalism as adopted in the literature for easy compari- 
son. First of all, there are 64 reactions among particles, 
whether atoms or molecules or free electrons. Each reac- 
tion is identified by a progressive number from 1 to 64 and 
so the associated reaction rate named fc„ (with n from 1 to 
64) . Each reaction will be a term in the righthand of the 28 
differential equations of system ([T]) and the associated re- 
action rate will coincide with one of the Rij terms of eqns. 
[1]). To establish the exact correspondence is a matter of 
patient work, which is not the prime i nterest here . Most 
of t hese reaction rates are taken from lAbel et al.l (|l997r) 
and Glover fc Savin to whom we refer for all the 

details. Some of these reactions are discussed in this sec- 
tion, and we pay attention to those rates that are vividly 
debated in literature. 

Minimal reaction network for CO. Since this 
molecule plays a key role in determining the properties 
of the ISM, we ha ve included a simp l ified network of re- 
actions taken from lNelson fc Langer (jl997t) to follow the 
creation/destruction of the CO in detail. These reactions 



of Woodall et al 



are listed in Table [5] tog ether with the corresponding rates 
(jioOTl ). Note that Tgoo = T/SOO, with T 



in K. 

Table [5] shows how CO is formed from C+ and O 
species, using CH, CH2, CHJ and CH;^ as intermediate 
products/reactants, and II2 and e~ as main reactants. 

Our model includes UV radiation, so we must consider 
the photodissociation of the molecules included in the net- 
work. The main effects of the radiation are on the carbon- 
based molecules and the corresponding photo-ionization 
rates are 

(4) 



rco = GnlO-i°s-i 



and 



r 



cm 



= Go5x 10-1" s^^ 



(5) 



where Go is the ratio of the adopted UV flux to that by 
Habing; i.e. /nab = 1-2 x 10^'* erg/cm^/s/sr. The expres- 
sion CHi^'' indicates that we are dealing with the CH, 
CH2, CH+, and CH+ molecules. 

Basing on eqn. ([S]), we must consider the products of 
all the ChI^^ molecules present in reactions like CHx'-^-' + 
7 — products. However, at this stage we are not inter- 
este d in the detail o f these reactions; therefore, follow- 
ing iNelso n fc Langeii (|l997t ). we introduce the parame- 
ter /3(Go) controlling the efficiency of the reactions from 
58 through 61 in Table [SJ This parameter is defined as 
l3{Go) = exp[Go • ln(0] with ^ = 5 x 10"!° from eqn. 
dSJ. In the absence of the UV radiation (Go =0), we get 
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Table 4. Reaction rates for the processes where the metals are involv ed a nd where cosmic rays (C R) in teract with 



Verner fc Ferlandl (Il99d) : (d) IVoronovl ()l997l) : (e) iGlover fc Jappsenl (j2007l) : (f) IZhao et al 



the IS M. References: . , . . ., . . , , 

(|2004[) : (g) IWalmslev et al.l (|2004[) . More details within the text 





Reaction 




Reaction Rate 
(cm /sj 


Rcf 


38 


C+ +e- - 


C + 7 


fcas = 'l>rcc(r, 6.85 X 10-^ 0, 11.3, 0.193, 0.25) 

A:39 = $coi(7;, 6.556 x 10""', 65.23, 2.446 x 10^ 0.7567) 


c 


39 


C + e" 


C+ + 2e- 


d 


40 


Si+ + e" 


Si + 7 


kio = $rec(r, 3.59 X 10"*, 0, 13.6, 0.073, 0.34) 

kii = $coi(Te, 8.616 X 10-1", 119.1, 4.352 x 10^ 0.7563) 


c 


41 


Si + e- - 


► Si+ + 2e- 


d 


42 


0+ +e- 


^0 + 7 


fc42 = $rec(r, 1.88 X 10-^ 1, 8.2, 0.376, 0.25) 





43 


0+e- ^ 


0+ + 2e~ 


fc43 = $coi(V,, 1.517 X 10~^ 360.1, 1.329 x 10^ 0.7574) 


d 


44 


Fe+ + e" 


-> Fe + 7 


k44 = $roc(r, 2.52 X 10-^ 0, 7.9, 0.701, 0.25) 


c 


45 


Fe + e" - 


¥ Fe+ + 2e" 


fc45 = 'I'coili;, 2.735 X 10-^ 1314.0, 4.659 x 10^,0.7568) 


d 



46 
47 



49 
50 
51 



0+ + H - 
+ H+ - 

+ He+ 

C + H+ - 
C+ + H - 
C + He+ 



0+ 

> 0- 



H+ 
+ H 

■ +He 



k46 = 4.99 X 10""T" 



fc47 = [1.08 X lO-^T'''^'^ + 4 X io-iOt'''°°'5''«] exp(-227/T) 



\ 0.3794 



-llTnO.517 

-^ ^ A J 

if(r < 6000 K) fe48 = 4.991 x 10"^^ (^j exp 

if(r > 6000 K) fc48 = 2.780 x 10"^^ ' ^ '^-^■^^^'^ 



( T \ 

\ 1.121x10" / 
T 



exp 



52 Si + H-* 



C+ +H 
C + H+ 
^ C+ + He 

> Si+ + H 



r 

14: I T 



53 Si + He+ Si+ + He 

54 Si + C+ Si+ + C 



fc49 = 3.9 X 10-16-^0.213 

fcso = 6.08 X 10-"* [iw""") ™P T- 

if(r < 200 K) fcsi = 8.58 X iq-i'^t"'^" 
if(200 <T< 2000 K) ksi = 3.25 x lO-i^T"'^*^® 
if(r > 2000 K) fcsi = 2.77 x IQ-^^T^ 
if(r < 10^ K) fc52 = 5.88 X 10" 
if(r > 10^ K) k52 = 1.45 X 10" 
fc53 = 3.3 X lO-'* 
ksi = 2.1 X lO-'' 



1^8.158x10^ y 



-19 jil.597 
T-13 jiO.848 
-13y 



55 
56 

57 



H + CR- 
H2 + CR 
H2 + CR 



H+ +e- 

> 2H 

> Hj + 



fcse = 0.46 CcR 
ksT = 1.50 CcR 
fc58 = 0.96 CcR 



Table 5. Reaction rates for the reactions belonging to the CO network. References: (a) IWoodall et al. ( 2007t ): (b) 
Petuchowski et al. I ()l989l) .See the text for more details. 





Reaction 




Reaction Rate 
(cm^/s) 


Rcf 


58. 


C+ + H2 - 


> CH++7 


fc58 = 4 X 10-1'^ Tgolj-" 


a 


59. 


CH+ + e- 


CH + H 


fc59 = 1.6 X lO-^Tgoj;" 


a 


60. 


CH+ + Ha 


^ ch;^ + H 


feeo = 1.6 X 10-* 


a 


61. 


CH+ + e- 


^ CHa + H 


fcei = 6.6 X 10-" 


a 


62. 


CH + ^ 


CO + H 


fc62 = 7.58 X 10-* T^oV 


a 


63. 


CHa + - 


> CO + 2H 


fc63 = 1.33 X W-^° 


a 


64. 


CO + c- - 


> C + O + e- 


/C64 = 2.86 X lO-^Tgog-^^ exp(-112700/r) 


b 



/3 = 1: this means that the existing photons cannot af- 
fect the formation of the various CHx^''. Instead, when 
Go = 1, it follows (3 — ^, which means that the forma- 
tion reactions have reduced their efficiency. Finally, when 
Go — ^ oo, (unphysical) then /3 — >■ 0, which means that 
the ChI'*''' molecules are completely destroyed before they 
can interact. The coefficient of the i-th reaction now is 
k[ = p ■ ki, with i e [58,61] and eqn. ([5]) is not included 
in the code. 

Photochemistry. In addition to this, we consider the 
group of photochemical processes listed in Table [5] for 
which we provide the cross section of the reaction and 



the analytical expression to derive the reaction rate as a 
function of the existing radiation field. 

In the present model of the ISM we do not include 
the chemical reactions with lithium and its compounds 



(e.g. L iH a nd LiII+) becau s e acc ording to iPrieto et al 



(|2008l ) and iMizusawa etal] (|2005[) . they are not impor 
tant coolants. After we neglecting lithium and compounds, 
the number of species to follow, including metals and deu- 
terium compounds, amounts to the 28 we have considered. 

There is some interest in the CO molecule because it 
is an important coolant at very low temperatures. Among 
others, recipes for the formation of the interstellar CO are 
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presented bv lRuffle et all (|2002l) and lGlover eTall (|2010t ). 
However, explicitly following the formation of the CO 
molecule would be too complicated for our aims. See below 
for the cooling rate by the CO molecule. 

In our chemical network we also include the effect of 
the cosmic ray radiation field. For convenience, we list 
in Table [H the i onizat ion rates for cosmic rays given by 
Walmslev et al. Cosmic rays are important for 

the cooling because they can destroy molecules that con- 
tribute to it, thus influencing the temperature of the in- 
terstellar medium (see the reactions shown in Table |4]). 
In particular, cosmic rays are able to destroy H2 and HD 
and to ionize atoms. The role of the cosmic rays is still 
largely unknown because the strength of their radiation 
field in various regions and epochs of the Universe is not 
firmly determined. While their effect can be explored in 
the vicinity of a strong source, for a random sample of the 
Universe there are no exact measurements of the cosmic 
ray radiation field. Furthermore, most NB-TSPH simula- 
tors like EvoL do not contain the full description of cosmic 
rays, but simply consider their presence as a parameter. 
Therefore, even if these reactions have been considered 
and their rates presented in Table |31 they have not been 
included in practice in the system of equations ([IJ. 

Photochemical reactions. For the photochem- 
ical reactions we ado pt the model and rates of 
Glover fc Jappsenl (|2007t ) and lVerner fc Ferlandl (|l996l ). In 
general the reaction rate is given by 



poo 

^photo = 47r / 



a{E)J{E) 
-E ' 



^(^)[l + /(i?)]di?, (6) 



where J{E) — J(hv) is the energy flux per unit frequency 
and solid angle of the impinging radiation field, cr{E) is 
the cross-section, t{E) the gas opacity at the energy E, 
and f{E) a numerical factor accounting for the effects of 
the secondary ionization, which are n egligible if the UV 
radiation is not dominated by X-rays ( Glover fc Jappsenl 
2OO7I) . The rate is in s'^. 

The cross sections <t{E) are expressed in two different 
ways according to the reaction under consideration. The 
reactions and associated rates are listed in TablelHl For the 
reactions from 1 through 8 we note that (i) the quantity 
e — \J {E/Et\i) — 1; (ii) for II2 and HD photodissqciation , 



we use the model proposed by [gI over fc Jappsenl (|2007l ). 
We have 



CTHa = 1-38 X l(fj{hv) (7) 

without considering self-shielding. For the HD we have 

CTHD = 1.5 X lO'^J{hv) , (8) 



X y 



0.5P-5.5 




(9) 



where y = \J •'^ f + ^f- ^^'^ meaning of the various sym- 
bols is given in Table [S) More details on the reactions 
and companion quantit ies an d symbols are g i ven b y 
Glover fc Jappsenl (|2007l) and IVerner fc Ferlandl (|l99d ). 
See also the references therein. 



Th e UV radiatio n field is calculated as in lEfstathiou 
1992h . IVedel et al.l (Il994h . and 



Navarro fc Steinmetz 



19971 ). In particular we have 



(10) 



where z is the redshift, i/jj = En/h is the frequency cor- 
responding to the hydrogen first-level energy threshold, 
ctTJV = 1 and J2i{z) is 



J2l{z) 



J2 



(11) 



l + [5/(l + ^)]4' 

where in our case J21 ~ 1. Finally, to derive the reaction 
rates, we integrate eqn. ^ from E = i?th to i? = cx) over 
the energy range of the radiation field photons. 

The integrals must be calculated at each time step if 
the radiation field changes with time or once for all at the 
beginning of the simulation if the radiation field remains 
constant. The integrals are calculated using Romberg's 
integration method. 

Heating by photodissqciation. Heating by pho- 
todissociation of molecular hydrogen and UV pumping, H 
and He photo-ionization, H2 formation in the gas and dust 
phase, an d finally ionization from c osmic rays are mod- 
eled as in iGlover fc Jappsenl (|2007l ). AU these processes 
are listed in Table [7] where is the critical density and 
Rd is the photodissociation rate. To describe the heating 
by the photoelectric effec t we adopt the model proposed 
bv lBakes fc Tielenj (|l994l) and discuss it separately in the 
Section 12.4.11 below. 

The H and He photo-ionizations increase the gas en- 
ergy at the rate 



r = 47r 



( 

E 



ME) (E-Eo)7j{E-Eo)dE, 

(12) 

where (t{E) is the reaction cross-section, J{E) is the pho- 
ton flux, e~'^(-^) is the optical depth for a photon of energy 
E, T] is the efficiency of the process (i.e. the fraction of en- 
ergy converted to heat), and {E — Eth) is the difference 
between the energy of the photon and the energy of the 
atomic ground level. The rates F are given in Table [71 



with hD = 12.87 eV. 

For the remaining reactions involving metals (i.e. 
from 9 through 12), th e rates contain the fits given by 
Verner fc Ferlandl (ll99(Th using the expression 



2.2. Reaction by reaction: notes on a few cases 

In this section we examine in some detail a few chemi- 
cal reactions that are widely debated and are affected by 
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Table 6. Cross sections of photochemical processes in cm^. The ene rgy E is in eV. The var i ous qu antities in use are 
desctribed within the text. References: (a) iGlover fc Jappsen ( 2007t ): (b) IVerner &: Ferland ( 1996 ). See these papers 
for further references. 





Reaction 
(cm''/s) 




Cross S6ctioii 


Note 


Ref 


1 

i 


-tl + 7 — ^ 


rl + e 


fji = D.o X iU ^ exp(4 — 4e — arctan s j [^i — exp — - — jj 


17^ 1 Q ft 
I^th = io.D 


a 


Z 


JJ + 7 — > 


JJ + e 




IT 1 Q ft 


a 


3 


He + 7 — 


¥ He+ + e" 


as = 3.1451 X 10"^^^^/^ [1.0 - 4.7416 ^^/^ + 14.82 ^ 


£th = 24.6 


a 








-30.8678 C^/^ + 37.3584 - 23.4585 ^^^^ + 5.9133 ^^j 






A 


n + 7 - 


rl -|- e 


-T. Oil V ^^^~^^^TT' t? ^3/2z7'— 3 

(T4 — - Z.ll X lU |^_c/ — -Lyth) ^ 


rjth — U. <oo 


a 


c 
D 


-112 + 7 — 


> rl + xl 


see text 






U 




-? rl -f- rl 


11 (^z.DO < H/ < ii.z(jfj5 — aex [— 4U.y ( + io.y/yoq 


IT n ftc: 

-C/th — Z.DO 


a 




















01 ^ fT'^OlA/T- /^(iv r Qfl OA -LVQCQCif^r 1 0001 A 

\i yvv.Li <_c/<zij (T5 — aex — ou.zD -f- / .oyoo ^ — i.zyzi4 












+0.065785 






/ 


-112 + 7 ^ 


> 1I2 + e 


11 (^10.4 < rl/ < iD.Oj fJei — y.ODU X lU ^ — y.4 X lU 


1 /I 

i^/th — io.4 


a 








if (16.5 <E < 17.7) fJ6 = 2.16 x 10"^^^ - 1.48 x 10"^"^ 












11 (ii./ <oU.UjfJ6 — i.ol X lU ^ — z./l 






Q 
O 




H + D 








9 


C + 7-^ 


C+ +e- 


0-7 = 5.027 X 10"^''$ph(1.607, 0.09157, a;/, 62.16, 5.101) 


Eth = 11.26 


a;b 










a;/(2.144, 1.133) 




10 


+ 7 ^ 


0+ + e- 


(78 = 1.745 X 10"^^$ph(0.1271, 0.07589, s/, 3.784, 17.64) 


£th = 13.62 


a;b 








(79 = 2.506 X 10-^^<l>ph (0.4207, 0.2837, x/, 20.57, 3.546) 


1/(1.240,8.698) 




11 


Si + 7 ^ 


Si+ + e" 


Eth = 8.152 


a;b 








(710 = 3.062 X 10~^^$ph(0.2481, 20.69, a;/, 2.671 x 10^,7.923) 


x/(23.17,rEsi) 




12 


Fe + 7 — 


Fe+ + e" 


EtH = 7.902 


b 










x/(0.05461, 138.2) 





Table 7. Heating processes. References: (a) iGlover fc Jappsen ( 200?!) : (b) IVerner fc Ferland ( 1996[ ). Details within 
the text. 



Process 



H2 photodissociation 
UV pumping of H2 
H photo-ionization 
He photo-ionization 

Gas-phase H2 formation 

Dust-phase H2 formation 
Cosmic-ray ionization 
Photoelectric effect 



Heating 

r = 6.4 X 10-'^iidnH2 



r = 2.7 X 10"" 
see text 
see text 



r = 
r = 



7.16 X 10-i2fcdustnnH (774^) 
r = 3.2 X 10-"Ctotn 
see Sect. 12.4.11 



Ref 



(2.93 X 10-i2fc2inH- +5.65 X V^-^'^kx^n^-^ (Tlfk;) 



large uncertainties. The reaction numbers are the same as 
those used in Tables [H O SI and E] for the sake of an easy 
identification. 

H2 + H+ -J> + H (reaction 9). Each author has 
his own favored rate for this reaction so that the overall 
uncertainty is large. Our pre scription is as follo ws. Up to 
3 X 10'' K we adopt the data of Savin et al. 1(12004'), based on 
the latest and most accurate quantum-mechanical calcu- 
lations of the vibrationally resolved cross sections for the 
charge transfer H2 + H+ -> + H at the center-of-mass 
collision energies from the threshold (^ 1.8 eV) up to 10 
eV. For higher collision energies, up to approximately 10^ 
eV, we derive the cross sectio n following t he suggestions 
by iBarnett et all (|l990l ) and Ijanev et all (,1987.) . which 



stand on the best kno wn experimental (see for instance 
Gealv fc van Zvll 1987 ) and theoretical data. , as we ll as on 
the recent measurements bv lKusakabe et all (l2003l). Th ese 
data were smoothly matched to those of Krstid (2002) at 
low energies, thus yielding the most updated cross section 
for the charge transfer from H to in the vibrationally 
ground state. This cr oss section i s show n in the top panel 
of Fig. m Following ISavin et all (|2004[) . from this cross 
section we calculate the rate coefficients for temperatures 
from the threshold up to 10^ K, as shown in the bottom 
panel of Fig. [T] 

We fit the reaction rate (in units of cm^s"^) with the 
analytical expression 
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4 5 
log[T(K)] 



Fig. 1. Top panel: charge transfer cross section for the re- 
action H2 -t-H+ +H. Bottom paneh Rate coefficient 
k{T) derived from the cross section in the top paneL The 
fits of k{T) (dashed hne) given by eqn. (0) and Table M 
cannot be visually distinguished from fc(T). 



log[A:9(r)] = ^aaog(r)\ 



(13) 



where T is the temperature in Kelvin. Two intervals are 
considered for the temperature, i.e. T = [10^, 10^] K and 
T = [10^ 10^] K and two different fits are derived. The 
coefficients for the two fits are given in Table \8\ The 
rates obtained from this analytical fit are compared to 
those derived from the numerical data. No difference can 
be noticed as shown in the bottom panel of Fig. [TJ 

Associative detachments of and D~ (reac- 
tions FROM 21 TO 24). These processes include the fol- 
lowing four reactions 

H- + H ^ H2 + e- , 
+H ^ HD + e" , 
+ D ^ HD + e" , 
+ D ^ D2 e~ . 

Following iGlover et al. I diooe!) and using their notation, 
the reaction rate is expressed as fc = const x ^ where 



Table 8. Fit coefficients for the charge transfer reac- 
tion H2 + H+ — > + H as described in the eqn. ([T^ . 
Coefficients are in the form a{b) = a x lO''. All the tem- 
perature ranges are shown. 





W <T < 10" K 


10" < T < lO** K 


ao 


-1.9153214(+2) 


-8.8755774(+3) 


ai 


4.0129114(+2) 


1.0081246(+4) 


a2 


-3.7446991(+2) 


-4.8606622(+3) 


aa 


1.9078410(+2) 


1.2889659(+3) 


04 


-5.7263467(+l) 


-2.0319575(+2) 


as 


1.0133210(+1) 


1.9057493(+1) 


as 


-9.8012853(-1) 


-9.8530668(-l) 


a? 


4.0023414(-2) 


2.1675387(-2) 



^ is a parameter taking values in the interval [0.65,5.0]. 
In our case we obtain fc2i — ^24 = 10~^ ^ and fc22 = 
^23 = 10~^^/2. In ROBO the parameter f can be varied 
in the above interval to investigate its effects on the overall 
results. In the present study we adopt ^ — 0.65. 

11+ + -> H -I- H (reaction 29). For the mutual 
neutraliz a.tion of H~ and H + we adopt the cross sectio n 
given by ICroft et al.l (Il999l) and iGlover fc Abel ^00^). 
The mutual neutralization rate is related to the rate of 
the associative detachment +11 ^ II2 +e^, in compe- 
tition for the available H" ions. B oth rates are impor- 
tant for the formation of II2 (see iGlover fc Abel I2OO8 . 
for more details). It is worth recalling here that other 
estimates for this re a ction rate have been reported by 
iGlover fc Abel! (|2008l ). lMoselev et al.l (|l970[) proposed the 
rate fc = 5.7 x io-6t-"-5 + 6.3 x 10"^ - 9.2 x 10-"T°-5 
-1-4.4 X 10^ ^"^T cm'^ s^^, whereas Dalearno fc Lepd (Il987h 
gave fc = 7.0 X IQ-'^T-O-^ cm^s'^ 



He^ 



He -t- 7 (reaction 35). This process 



can occur either via direct radiative recombination or di- 
electronic recombination, followed by radiative relaxation. 
Therefore, the reaction rate is the sum of two terms 

fc35 = fc^ + fc§5 , (14) 

in units of cm'^s"^. The first term is the di-electronic re- 
combination rate, described by 

fc^5 = 1.544 X 10-9r-^-^exp(-48.596/Te) 

X [0.3 + exp(8.1/re)] , (15) 



where is the temperature expressed in eV. The second 
term is the radiative recombination rate whose tempera- 
ture dependence is 



fc^5 = 3.925 X IO^^^T; 



0.6353 



(16) 



Both terms are taken from lAbel et al.l ([1997|). 

Z+ + c" ^ Z + 7 (reactions 38, 40, 42, 44). An 
important process to include is the metal recombination. 
The metals considered by our ISM model and ROBO are 
C, Si, O, and Fe. The recombination rates of the metals 
are calculated ac c ordin g to the formalism proposed by 
Verner fc Ferlandl (|l996l ) 
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^rcciT,A,To,Ti,b) = 



To(l 



To)'"' (1 + ^1)'+'' 



(17) 



in units of cm'^s ^, where T is the temperature, Tq = 
V'T/to and Ti = ^Tjr^. In Table g] we give the fit coef- 
ficients A, 6, To and t\ for each metal. 

Z-he" ^ Z+-h2e- (reactions 39, 41, 43, 45). For 
the collision s between electr ons and metals we use the fit 
proposed bv IVoronovl ( 1997 ) expressed by 



$coi(T,A,P,A^,X,if) = A 



(18) 



in units of cm'^s^^, where I] = /S.E/T in which /S.E is the 
energy difference between the two atomic levels involved 
in the process, and T is the temperature. Both Ai? and T 
are in eV. The parameters P, AE, X, and K are given 
in Table H 

2.3. The dust 

Dust grains take part to the process of molecule formation, 
e.g. HD and H2 form on the surface of dust grains; there- 
fore, all the physical processes involving the dust must be 
described and treated with the highest accuracy. To un- 
derstand how grains take part in the process of molecule 
formation, we need to know the mechanisms of grain for- 
mation and destruction, along with the distribution of the 
grain temperature and size. 

Given an initial set of dust composition and abun- 
dances, our dust model follows the evolution of the dusty 
components during the history of the ISM due to the cre- 
ation of new grains and the destruction of the existing 
ones. Creation of new dust grains is governed by several 
processes that have different efficiencies depending on the 
size of the dust particles. The same applies for the destruc- 
tion that is mainly due to the thermal motion of the gas 
particles and shocks from supernovas. Thermal destruction 
is quite easy to model, because the only parameter at work 
is the gas temperature. Shock disruption is more difficult 
to evaluate. The main uncertainty comes from discrete 
numerical hydrodynamical simulations only being able to 
follow shocks up to a given (often too coarse) resolution, 
yet insufficient for the microscopic description required 
here. To cope with this, we have followed a "statistical" 
approach. In brief, once identified the gas particles of the 
NB-TSPH simulations with their turbulent velocities (sug- 
gested by their velocity dispersion), we assume that all the 
shocks inside them follow the Kolmogorov law. See below 
for more details. 

Furthermore, grain destruction may depend on their 
size. This is the case of the destruction by thermal mo- 
tions, where the small grains are destroyed before the large 
ones. As a consequence of this, the size distribution func- 
tion and the abundances of dust grains of different type 
are continuously changing with time. The ever-changing 



size distribution function plays an important role in the 
formation mechanism of HD and H2, which are among 
the most efficient coolants. In brief, changing the distri- 
bution function of the dust grains means changing the 
q uantity of key-role molecu les produced by dust, as shown 
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We analyze here the different aspects of the grain evo- 
lution and their role in the formation of coolant molecules. 
First, we focus on the distribution function of the dust 
grains, then we describe the formation of coolants on 
grains. Finally, we analyze the destruction and the forma- 
tion of dust, as well as the effects of the grain temperature. 

2.3.1. Size distribution function of dust grains 

We adopt a simp le power-law, MRN - like, grain size distr i- 
bution function (jMathis et al.lll977t iDraine fc Ledll984l ). 
This is given by dn(a)/da oc a"'*', with a the grain di- 
mension, n(a) the corresponding number density of grains 
with dimension a, and A — 3.5; this distribution covers 
the range 5 A < a < 2500 A and is extended to the 
smallest grain dimensions in such a way as to include the 
typical polycyclic ar omatic hydrocarbons (PAHs) region 
( Li fc Drainell2001bl ). Eve n if more complicated distribu - 
tions have been proposed ( Weingartner fc Drainel[2001a ). 
our simple choice is suitable as an initial condition for the 
purposes of this work. 

Indeed, it is worth pointing out here that the initial 
size distribution function changes with time owing to dust 
destruction and formation, which in turn depend on the 
grain size. Consequently, the power law we started with 
does no longer apply. This will soon affect the formation 
of H2 or HD. For this reason, the knowledge of the current 
size distribution function is paramount. 

First of all, we split the interstellar dust in two main 
components: the carbon-based (thereinafter simply car- 
bonaceous grains) and the silicon-based grains (there- 
inafter the silicates). These play different roles in the for- 
mation of molecular hydrogen, as well as in different for- 
mation and destruction rates. In principle, there should 
be an additional group to consider, namely the PAHs, be- 
cause the three types of dust have different efficiencies in 
the H2 and HD formation. However, there is the lucky 
circumstance that PAHs and graphite grains have similar 
efficie ncies to those shown by Fig. 2 of Cazaux fc SpaansI 
Basing on this, we lump together graphite grains 
and PAHs and treat separately the silicates. Therefore, 
thereinafter we refer to the mixture graphite grains plus 
PAHs as "carbonaceous grains" . 



2.3.2. Dust-driven H2 and HD formation 

The model wc adopt t o describe the so-called d ust phase 
is the one proposed by Cazaux fc SpaansI ( 2009t ). In brief, 
the dust phase has its own network of reactions that estab- 
lish the relative abundances of the various types of dust 
and govern the formation of H2 and HD. We can neglect 
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the presence of other molecules like CO. The newly formed 
H2 and HD are fed to the gas phase and enter the system 
of equations governing the abundances of the ISM and 
vice versa. In reality, the two phases should be treated si- 
multaneously, governed by a unique system of differential 
equations determining the number densities of the elemen- 
tal species, dusty compounds, and free electrons at once. 
In practice this is hard to do because in general the time 
scales of the various processes creating and destroying the 
elemental species of the gas phase are much different from 
those governing the formation / destruction of dust grains . 
To cope with this difficulty, the ICazaux fc Soaand (I2OO9I) 
model separates the two phases, provides the results for 
the dust phase, and more important provides a link be- 
tween the two phases that allows determining the for- 
mation rate of II2 and HD by dust as a function of the 
gas temperature Tg and dust temperature T^. This link 
is named the "Sticking Function S{Tg,Td)" ■ The sticking 
function somehow quantifies the probability that an hy- 
drogen atom striking a grain sticks to the surface rather 
than simply bouncing off. The sticking function charac- 
terizes the probability for an atom to remain attached to 
a grain. Here we stri ctly follow this way of pr oceeding. 

According to the lCazaux &: SpaansI ( 20091 ) model, the 
rate at which H2 molecule forms over the grain surface is 



i?d(H2) = -n(H)wHndaeH.5(Tg,Td) cm-^s"!. (19) 

where n(H) is the density of hydrogen in the gas phase, 
Vh = \J%kTg/ (Trmn) is the gas thermal speed (fc is the 
Boltzmann constant and mu is the hydrogen mass), the 
dust number density, a the grain cross section (i.e. tto^), 
the intrinsic efficiency of the process, and S{Tg,Td) 
the sticking function. This latter is in turn a function of 
the dust and grain temperatures 



S{Tg,Td) = 



H-0.4 



100 



H -1 



0.5 



0.2|i| (20) 



where Tg and Td are in Kelvin. With this function the 
probability for a gas molecule to remain stuck on a dust 
grain is higher in a cold gas than in a hot one. This prob- 
ability is even higher for the cold grains. For HD the ex- 
pressions are similar. 

Equation (jl9p contains the intrinsic efficiency (or 
as eHD for HD), which is the probability for the process to 
occur. In our case it is the probability that atoms, which 
are stuck to a grain, react to form H2 or HD. 
For the carbonaceous grains, the efficiencies ejja and epjo 
coincide and are equal to 



+ 0-25(1 + /^^) 



(21) 



Table 9. The values used to compute the formation effi- 
ciency of H2 and HD, with j?phv, Ec\i , and Eg in K, Cpc 
is in A. From Cazaux fc SDaansr ( 2009h . 



Surface 


-Ephy 




Es 




Carbons 


800 K 


7000 K 


200 K 


3 A 


Silicates 


700 K 


15000 K 


-1000 K 


1.7 A 



The efficiency for the silicates esi is 



esi 



1 



16Td 



ape \/ -Bphy— -Es 



-F, (22) 



where /3d = 4 x 10^ for H2 and fid = 5.6 x 10^ for HD. 
The term is a function of the gas temperature and can 
be written as 



F{T) = 2 



(23) 



where Th is given by the expression 



Th = 4 




Ephy — Es 



exp 



E, 



ch 



Es 



E, 



ch 



(24) 

The various quantities appearing in the above relation- 
ships are listed in Table El For mo re details on these equa- 
tions see Cazaux fc Soaansl ( 20091 ). Comparing eqns. (^11) 
and (j22p . we see that the efficiency is high when the dust 
temperature is low. For the silicates the efficiency window 
is shorter than for the carbonaceous grains. For the sili- 
cates the efficiency is high for temperatures up to 20 K 
and then falls by two orders of magnitude. The carbona- 
ceous grains are efficient for temperatures up to 100 K, 
where the efficiency is still 0.1 (instead of 0.01 as for the 
silicates). Finally, it is worth noticing that the efficiency 
profile is smoother for the carbonaceous grains. 

From all these considerations it follows that cold dust 
and warm carbon-dominated dust in a cold gaseous envi- 
ronment are the most efficient drivers for the formation of 
coolant molecules like H2 and HD. 

2.3.3. Grain formation 

Here we briefly examine the formation of dust grains. We 
start with an initial number density with the size distribu- 
tion given in Section [2.3.11 and with a given ratio between 
the silicates and the carbonaceous grains. The latter is 
a free parameter varying in the range [0,1], where zero 
stands for dust made of sole carbonaceous grains; one is 
for dust made of sole si l icates . 

According to iDwek (Il998l) . the temporal variation in 
the size distribution function of grains caused by accretion 
is given by 
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dn{a) 
dt 



Cd a{Tg, Td)TTa^ngnd{a)vg 



(25) 



where all the quantities have the same meaning as in eqn. 
(frail, but for a{Tg,Td) and c^. 

The quantity a{Tg,Td) is a sort of sticking coefficient 
depending on the gas and dust temperature and the type 
of dust. This coefhcient is therefore forced to change in the 
course of the evolution. Furthermore, (and consequen- 
tially n) are functions of a. Equation (j25l) is similar to eqn. 
(fT9)) because the process is similar, except that now this 
mathematical description is applied to the carbon atoms 
that remain stuck to the carbon lattice of the grai n. Our 
expres sion differs slightly from the original one of iDwek 
( 1998h because the parameter Cd is in troduced in e qn. 
to take some considerations made by Dwek (JlOOS") himself 
into account. Inserting Cd = 1 we can obtain the time 
scale of the process. In a standard cold gas this time scale 
is r 2 X 10^ yrs, which is significantly smaller than the 
observational estimates. Normal evaporation, caused by 
cosmic rays or UV heating and grain-grain collisions can 
halt the growing of the dust grains. The factor Cd somehow 
takes this phenomena into account. It is estimated to be 
on the order of Cd = 10^'^. Wc name Cd the delay factor. 

For the sticking coefficient a{To,Td), we make use of 
the data bv lLeitch-Devlin "fc WiUiamal (|l98,5h and consider 
a carbon atom of the gas phase as impinging on a carbon 
lattice. The equation fitting the data has the form 



a{Tg,Td) = 0.0190 (0.0017 Trf 
X exp (-0.0070 Tg) , 



' 0.4000) 




1000 



Fig. 2. The 

molecules 



sticking coefficient for the carbon 

over a ca. rb on lattice derived from 

Xeitch-Devlin fc Williamsl (|l985[) . The different lines 
indicate a different dust temperature: 3 K (solid), 100 K 
(dotted), and 200 K (dashed). Tg is the gas temperature 
in Kelvin. 



(26) 



stars, and remnants of supernova (jDweklll998[ ). The rea- 
son behind this is that all of these are external sources 
of dust that eventually enrich the ISM in dust content 
and therefore determine a different initial dust content as 
input for ROBO. Indeed, for our purposes, only the in- 
ternal sources of dust in the unit volume are important. 
The information about the initial conditions of the dust 
mixture should be provided to ROBO by the companion 
NB-TSPH code EvoL, of properly taking the dusty yields 
coming from the stars into account. 



where Tg and Td are the gas and dust temperatures, re- 
spectively (see Fig. [2]). The sti cking coefficient a{To,Td) 
has no dimensions. The data of iLeitch-Devlin fc Williams 
(19^ ])rovide a good dependence on the gas temperature 



Tg, but a poor one on the dust temperature Td- Therefore 
an accurate fit is not possible. The above relationship can 
be used in the temperature intervals 10 K < < 1000 K 
and 3 K < Trf < 300 K. At present we also use the same 
model for silicate grains even though this approximation 
might not be accurate. We plan to improve upon this point 
in the future. 

With this model the formation efficiency is higher 
when the interstellar medium has a temperature T 
100 K and when the dust grains have a temperature of 
about 300 K. For dust grains with temperature higher 
than 300 K or lower than 3 K, the formation efficiency is 
unkno wn. Though the fits from lLeitch-Devlin fc Williams 
(|l985h are not very accurate, they are accurate enough for 
our purposes because we are only interested in the shape 
and maxima of the curves in Fig. [2l 

Finally, we note that we have already described the 
formation of dust by means of a general process of accre- 
tion in the ISM, leaving aside other sites of dust forma- 
tion like the envelopes of obscured AGB stars, Wolf-Rayet 



2.3.4. Grain destruction 

In our model we assume that there are three processes 
destroying the grains of the ISM, i.e. shocks, in particular 
supernovae induced shocks, vaporization by high-velocity 
shocks, and thermal sputtering, in which dust is destroyed 
by the thermal motions of the gas. We need to describe 
the three processes in a way suited to the NB-TSPH for- 
malism. 

Destruction by shocks. This phenomenon is diffi- 
cult to model. First of all shocks may easily induce turbu- 
lence in the ISM and to be suitably described one needs 
some assumptions about the velocity fields of the fluid 
elements. The main difficulty arises from the fact that 
within an SPH gas particle, because of the typical mass 
resolution, unresolved shocks may take place at different 
velocities. To avoid this difficulty, we treat the g as as a 
turbulent fluid, assuming that even inside a gas parti- 
cle the turbulent nature of the fluid is preserved, and use 
the Gaussian velocity distribution (f>{v) derived from the 
Kolmogorov law for the power spectrum E(k) oc k~°'dk 
with a=5/3. Second, grains of different masses and sizes 
move at different velocities: the large, more massive grains 
moving slower than the smaller less massive ones. This 
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suggests adopting the point of view in which the small 
fast grains are considered as projectiles impinging on the 
massive slow grains considered as targets. To summarize, 
the ISM we are dealing with is turbulent because of the 
shocks crossing it and rich in dust grains of all possible 
dimensions and masses interacting with each other and 
with the shock fronts. 

Given these premises, we model the destruction of 
dust grains by shocks according to the picture by 
Hirashita fc Yan The dust grains are grouped in 

N = 100 bins of number density according to their mass. 
The mass of a grain is in turn expressed as the product 
of a reference density, here assumed to be the density of 
the graphite grains pgi- = 2.3g/cm'^ times the volume of 
the grain Vgr — 4/37ra^ where a is the radius of the grain 
(assumed to be spherical for simplicity). Therefore each 
density bin contains grains whose mass goes from miow 
to rUup, where both mass limits vary with the bin. The 
corresponding radii are given by miow — 4/37rafQ^ and 
— 4/37ra^p. The grain radii follow the distribution 
law given by dn{a) = Cnor x a~^-^da, where Cnor is a 
suitable normalization factor to be determined. Therefore, 
the number density of the grains in each mass interval is 



n{a) da ~ Cno 



-3.5 



da 



(27) 



The total number density of dust grains of any size and 
hence mass is 



1=0 



(28) 



When a projectile of mass nij hits a target of mass to^, 
the target loses a fraction of mass /shWi if TUj < mi/2, 
or it is totally destroyed if nij > mi/2 (i.e. fsh = !)• In 
both cases the projectile is destroyed. The lost mass, i.e. a 
fragment of lower mass, so the dimension is assigned to the 
mass (size) bin according to the adopted power law. The 
remaining part of the target (remnant) of mass {1 — fsh)mi 
is added to the appropriate mass, hence size bin, so the 
number density variation of the i-th bin is 

where the subscript L indicates the term representing the 
mass lost in fragments, F the mass gained from the frag- 
mentation, and finally R the mass of the remnants moving 
from other bins to the i-th one. The first term must be 
negative. 

To estimate the various contributions to the total num- 
ber density variation, we must consider the probability of 
an impact. This will be proportional to the number den- 
sity of the targets n^, the density of the projectiles rij, 
their relative speed u^, and their sizes (a^ and aj). The 
assumptions on the speed will be discussed later in this 
section. For the targets belonging to the i-th bin, we find 



that the mass lost in fragments corresponds to the density 
change of 

j<i 



dt 



^ij ^ij^i^j , 



where 



Ui+rij 
1 



if mj < mi/2 , 
if mj > mi/2 



and the impact coefficient a^jis 

a^j = VijTT{ai + Oj)^ , 



(30) 



(31) 



(32) 



Consequently, the total density variation of newly created 
fragments is 

2^ r\i ■ 



dt 



(33) 



They will be distributed among the bins with a power law 
to obtain the (dni/dt)p of each bin. 

The variation in the remnant number density is given 

by 

'dnk\ , .dn, 

R 



dt 



dt 



where k is the index of the bin receiving the remnant given 
by 

mi(l - fsh) - TOio- 



k = int 



Am 



(35) 



with i the index of the initial grain that suffered frag- 
mentation and Am the mass range of each bin, {m^p ~ 
m\o-w)/N. 

The shock velocities obey a Gaussian distribution (w) 
normalized to 



4>{v) dv 



1 , 



(36) 



where wiow and Wup are the limits of the shock velocities 
over which the normalization of 0(w) applies, thus fixing 
the normalization constant. The relative speed Vij of the 
grains is obtained from the velocity distribution of eqn. 
pop with flow = 1 km/s and = 200 km/s and a total of 
30 velocity bins. The velocity 0(w) is a Gaussian centered 
&tv = 100 km/s with a dispersion of 30 km/s. The integral 



4>{v)dv , 



(37) 



yields the relative weight of each velocity bin. The relative 
velocities of the projectiles are distributed according to 
these weights and the total density of targets in the i-th 
bin changes as described by eqns. ([50]) and ([5^. 

Vaporization. Since we also deal with high-velocity 
shocks, vaporization of dust grains into the gas phase may 
become important. We use the formalism of Tielens et al.l 
(jl994al ) for the impact of carbon grains. The fraction of 
vaporized material is 



vap 



(38) 
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with fyap given by 



/v. 



fvl + fv2\r^ 



(39) 



where vt = 23 km/s is the shock threshold velocity, nii 
and rrij are the masses of target and projectile grains, 
respectively. Equation ([55]) must satisfy the condition v > 
Vt- The quantities fyi and /„2 are taken from lTielens et aL 
(Il994al) and are given by the relations 



fvl = 2.59- 



and 



where 



fv2 



2.11 

8/9 



Tl 



0.3 (s + X-i-0.11 



1.3 



(40) 



(41) 



(42) 



(s + X-i-1) 

with M the Mach's number and s — 1.9. 

Destruction by thermal sputtering. To evaluate 
the fraction of grains destro yed by thermal sput t ering , we 
adopt the approximation by iDraine fc Salpeter ( 1979t ) . A 
grain of dust in a medium with temperature T < 10^ K 
has a destruction time (i.e. its lifetime) 



a certain mass function, hence mass-luminosity law) are 
born. Without coolants neither structures nor stars would 
form. In this context, the grains and their temperature 
in turn play the key role in the formation of efhcient 
coolant molecules like H2 and HD. Cold grains form more 
molecules than the warm ones. Since the dust tempera- 
ture depends on the surrounding photon flux, it means 
that stars in the neighborhood are crucial for heating the 
dust particles. All this forms a closed loop of interwoven 
physical processes: dust forms coolants - coolants induce 
star formation - stars heat the dust. Understanding the 
details of this mutual interaction will allow us to get clues 
on the star formation history in general. 

In this scene a starring actor is the gas cooling. 
Chemical reactions are sensitive to the gas temperature, 
hence to the gas cooling; indeed, the formation of coolants 
depends on the gas temperature, which in turn depends 
on the cooling process. We split the cooling of the gas in 
several sources characterized by the physics of the dom- 
inant process. Above lO^'K, there are two very p opu- 
lar descriptions of the cooh ng, namely |Ce3 (119921) and 
Sutherland fc Dopital ( 1993). At lower tem peratures we 



have the metal cooli ng (iMaio et all 2007 ) , the molecu 
lar h ydrogen cooling ( Glover fc Abe l 2008; Galli fc Palla 
19981 ). and finally, the HP cooling (|Lipovka et al.ll2005l) . 



Tdist ~ lO^a/nn yr. 

with a in nm. From this we can obtain the destruction 
rate per second. The above relation is only valid for high 
temperatures. 

To include the ter nperature depe i idence of the destruc- 
tion rate, we refer to iTielens et al.l (|l994bl ). with the aid 
of which we model the dependence on the gas temperature 
of the lifetime of the dust. As expected, dust grains with 
lower temperature have a longer lifetime. Finally, accord- 
ing to Draine fc Salpeter ( 1979f) . the small grains have a 
shorter lifetime than the large ones in the same environ- 
ment. 

The grain temperature depends on the radiation field 
generated by all the stellar sources. However, for the sake 
of simplicity, in this study we used a fixed value for the 
grain temperature. Therefore, a big improvement would 
be gi ven by including a d escription of the photon diffusion 



(e.g. iMathis et al.lll983l) . Our choice is partially justified 



by the fact that the companion NB-TSPH code EvoL 
does not yet include photon diffusion. This leads us to 
postpone the implementation of different grain tempera- 
tures to when photon diffusion is included in EvoL. 

2.4. Heating and cooling 

At this stage, it is worth discussing in some detail the 
role played by the heating and cooling during the history 
of star formation in a galaxy or a cosmological simula- 
tion. We have already emphasized that coolants are the 
key elements for the formation scenario. As the gas cools 
down, more and more spatial structures and stars (with 



2.4.1. Heating by photoelectric ejection of electrons 
from dust grains 

The photoelectric ejection of electrons by dust grains 
is an important source of heating. The model pro- 
posed here is based on iBakes fc Tielend (|l994[) and 
IWeingartner fc Draine ( 2001b ) . The photoelectric heating 
is given by 



(43) 



H{Nc.Z) = WtI / aab.(iVc)llon(^C,/Pz) 



where W is the FUV dilution factor, Uabs the photon ab- 
sorption cross section, l^on the photoelectric ionization 
yield, F(v) the UV radiation flux, and g{NcTlPz) the 
kinetic energy partition function. Here, Nc is the number 
of carbon atoms that form the dust [Nc ^ 0.5 a'^ with 
a the r adius of the grain in A, according to Li fc Draind 
(|2001bf )) and IPz the ionization potential of a grain of 
charge Z, and vh is the Lyman frequency and vz the fre- 
quency corresponding to IPz- In our model the range of 
grain sizes goes from 5A to lOOA , since the total heatin g 
is mostly due to small grains (see iBakes fc Tielens 1994 , 
for details). 

The total photoelectric heating is 

r= / y^H{Nc,Z)f{Nc,Z)n{Nc)ANc , (44) 
Jn^ z 

where f{Nc, Z) is the probability to find a grain composed 
of Nc carbon atoms at a certain charge Z. This can be 
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computed by considering the collisions with electrons and 
ions. The detailed balance yields 

fiZ) [Jpe(Z) + Jion(^)] = fiZ + l)Je{Z + 1) , (45) 

with Jpo the rate of photoelectron emission, Jiom and 
Je the accretion rates of ions and electrons. For the de- 
tailed calculation of f{Z) see 
while for all the details ab out 
Weingartner fc Draind (|2001bl ). 



Bakes fc Tielens 



Jue^ J'lon and 



(|1994D . 

Je see 



As the FUV absorption cross section dabs depends on 
the absorption e fficiency of FUV phot ons b y the grains, we 
use th e data of Draine Sz Lee ( 1984[) and Laor fc Draine 
(|1993[) for the graphite grains andl Li fc Draind (|2001al lby 
for the PAHs. We also use the model proposed by 
Li fc Draine (I2001allbl l3) to create a mixed population of 
PAHs and graphite grains. Following their paper the total 
absorption efficiency is 



Qabs(a, A) 
defining 



?rsKA), (46) 



f (a) = (1 - q)miii 



(47) 



by 9 = 0.01 and aq = 50 1 fsee lLi fc Draindl2001allbl l3. for 
more details). 

In the same way we can also obtain the cooling asso- 
ciated with the electron recombination. This is given by 



C{Nc,Z)^n,s, 



8kT 



Trnii 



7ra^A(T, v) , 



(48) 



where rii is the density of the charged particle with sticking 
coefficient Si — 1 and mass rui and a is the grain size. 
The te rm A(t, v) is described in detail in lBakes fc Tielens 
(Il994l) . The total cooling is obtained as in eqn. (H^ . In 
Fig. |3] we show the results from a simple model with a 
grain distribution with ridust = 10^^ cm~^ and a gas with 
a temperature of 10^ K. 



2.4.2. Which are the sources of cooling to consider? 

Cooling at high temperatures (T > 10*). Two main 
sources of cooling can be found in th e literature: (i) The 
formulation proposed by CenI ( 1992h which includes the 
following mechanisms: 



— collisional ionization - H, He, He^, He(2'^S), 

— recombination - H+ , He^ , Hc^^ , 

— dielectronic recombination - He, 

— collisional excitation - H(all n), lle^{n ~ 2), He (71 ~ 
2,3,4), 

— bremsstrahlung - all ions, 

which is particularly suited in presence of non equilibrium 
re actions for the H and H e chemistry, (ii) The description 
of [Sutherland fc Dopital fl993.. thereinafter SD93), which 
is particularly suited to treat cooling in presence of metals 
and to describe it as function of the gas metallicity. Their 
case for Z = can be left aside. Both cooling models 
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Fig. 3. Top panel: heating by dust (solid line), cooling 
(dotted line) and the difference between heating and cool- 
ing (dashed line) for different diameters of the dust parti- 
cles. Bottom panel: the probability of finding a grain with 
a certain charge Z as a function of the grain dimensions. 



are used for temperature higher than 10'* K. We refer to 
them by the acronyms C92 or SD93. In our models of the 
ISM we adopt both sources as appropriated to the current 
physical conditions. 

A difference between C92 and SD93 is that the cooling 
rate is described in the former by numerical fits (no further 
interpolations are needed), whereas in the latter one has 
to interpolate huge tabulations of data in temperature and 
metallicity. To this aim we used a surface fit (routine SFIT 
in IDL) where the plane surface passes through four points 
that correspond to two discrete values in temperature and 
two discrete values in metallicity. The analytical expres- 
sion of this plane yields the cooling rate. Given the two 
pairs of interpolating points, {Tq,Zq), {Tq,Zi), {Ti,Zq), 
and (TijZi), the analytical function that describes the 
surface can be written in the form 



A; 



SD93 



{T,Z) =ao + aiZ + a2T + a3TZ, (49) 



with tti the fit coefficients. So, for a generic point (T, Z), 
for which Tq < T < Ti and Zq < Z < Zi, the rate cooling 
is given by eqn. (|49l) . 
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H2 co oling. The H2 cooling requ ires a different de- 
scription. iHollenbach fc McKed (|l979l ) first evaluated the 
molecular hydrogen cooling. A mod e rn, w idely used source 
of H2 cooling is by iGalh fc Pallal (|l998l ). Although this 
cooling function is quite accurate, we prefer to add 
some more de t ails b y including the function found by 
Glover fc Abel (|2008l ). which takes not only the H2 - H 
interaction into account, but also the collisions with He, 
H+, H2, and free electrons. It is described by 

Ah2 = X! X! ^'^^ t"**^ logics)'] nn^Uk ergcm"^ s~^ , 

k i 

(50) 

where T3 = K with T the gas temperature, 

aik is the i-th fit coefhcient of the k-th species {k = 
{H, H+, H2, e~. He}), and n are the number densities. The 
orto-para ratio is assumed to be 3 : 1. Outside the tem- 
perature r ange of the fit s we u se the molecular hydrogen 
cooling bv lGalh fc Pallal (Il998h . 

HD cooling. To describe the cooling by the deuter- 
at ed molecular hydrog en we use the model proposed 
bv iLipovka et al. (j2005[l . It includes HD roto- vibrational 
structures, radiative and coUisional transitions for J < 8 
rotational levels, and the vibrational levels w = 0, 1, 2, 3. It 
has been found that including the roto-vibrational tran- 
sitions increases the cooling efficiency of the HD. The fit 
provided by the authors depends on the gas density and 
temperature. It can be parameterized as 



A(T, n 



HD 



EE' 

1=0 j=0 



logio(n)Mogio(T)^ 



(51) 



where Cy is the matrix whose elements are given in Table 



1101 Using eqn. (1511) is particularly convenient from a nu- 
merical point of view as it provides fast evaluations of the 
cooling by the HD molecule. 

Cooling by the CO molecule. Owing to founding 
hypotheses of our model of the ISM presented in the in- 
troduction, the ISM is optically thin so that it is worth 
considering the straight cooling of the CO molecules by 
radiative transitions among rotational and vibrational en- 
ergy levels. 

To t his aim, we may a dopt the steady-state equa- 



tions of McKee et al.l ( 1982h and the rates calculated by 
Schinke et al.~(ir985l) by appl ying the L coefficients 
to the iGoldflam et al. 1 (Il977h method and including the 
Depristo et al. I (I1979I) correction (see these studies for 
more details). 

To calculate the total rotational coo ling, we use 



the eq uation of statistical equilibrium from iMcKee et al, 
(1983): 



TlH: 



E 



2 / , 



0, 



(52) 



where rij is the number density of the CO molecules in 
the state j — J, the number density of the molecular 



hydrogen, 7^ the rate coefficient of the transition, and Aj 
the A-coefficient from the state j to j — 1. The system 
of equations is completed by the condition J2j "-j — "-co • 
The solution of this linear system is straightforward. The 
emission due to the transition from J to J — 1 is 

Ij — h^L'^jI—L erg/s/molecule/sr, (53) 
47r nco 

from which we get the total cooling rate Aiot — 
nco J2 jlj- 

The vibrational component is taken from 
IHollenbach fc McKed (|l979l) . considering collisions 
with H2 whose rates, are 



701 = 4.3 X 10"i''Texp(-68/T3)exp(-Sio/fcr) 

702 = 2.7 x 10-i3Texp(-172/T3)exp(-£;2o/fcT) , (54) 

where Eio/k = 0.5£'2o/fc = 3080 K, T3 = T^/^ and k is 
the Boltzmann's constant. 

Cooling by the metals. The cooling process by the 
metals is included using the list of the metals adopted by 
iMaio et al.l (l2007l). namely CII Sill, 01 and Fell (see the 
Appendix B in lMaio et al.ll2007() . This source of cooling is 
particularly important for temperatures lower than 10^ K. 
The cooling is due to the fine structure level transitions 
of the ionized carbon 2P(3/2) — >■ 2P(l/2), and similarly 
for the ionized silicates. There are two other species that 
take part in this process: neutral oxygen (nine transitions 
between levels ^Sq, ^D2, ^Po, '^Pi, and ^P2) and ionized 
iron (five transitions between levels ^D^ with i odd in the 
range [1,9] G N). The cooling from each transition is ad- 
ditive and is given in erg/cm^/s by 



A. 



jfi -|-7|i«'e/nH 



7" + 75 + (7-, + 7° ) nc/nn + A,,/n^ 



(55) 



where jj^ is the reaction rate for the hydrogen, and j°j is 
the same but for the electrons, Aij is the Einstein's coef- 
ficient between the i-th and j-th levels, AEij is the energy 
difference between the levels, ntot is the total number den- 
sity of the species considered, and finally, n^i and rie are 
the hydrogen and electron number density, respectively. 
To complete the equation we need to know 



7i^ 



9j 



exp [-AEij/ (fcsTgas)] 



(56) 



both for hydrogen and electrons; gi and gj are the level 
statistical weights, ks is the Boltzmann's constant, and 
Tgas is the gas temperature. 

If the model is calculated at high redshift, the CMB 
pumping must be considered. This process is described by 
including the stimulated emission coefficient in the treat- 
ment of the cooling. The black-body radiation field of the 
CMB has a temperature that depends on the redshift as 
Tcmb{z) = (1 + z)Tcmb{0) K- However, as the models of 
the ISM refer to the current age, the temperature is kept 
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i = 


i = 1 


i = 2 


i = 3 


i = 4 


j 


= 


-42.56788 


0.92433 


0.54962 


-0.07676 


0.00275 


3 


= 1 


21.93385 


0.77952 


-1.06447 


0.11864 


-0.00366 


3 


= 2 


-10.19097 


-0.54263 


0.62343 


-0.07366 


0.002514 


3 


= 3 


2.19906 


0.11711 


-0.13768 


0.01759 


-0.00066631 


3 


= 4 


-0.17334 


-0.00835 


0.0106 


-0.001482 


0.000061926 



constant and equal to the present day value (see below). 
In any case, owing to the low value of the present-day 
CMB temperature, this effect is neglected. To obtain the 
total cooling by the metals we use Amctais = X^i where 
i indicates a generic level transition for any of the four 

metals we have included. 

The cooling proposed by iMaio et al.l (j2007l ) 
is improved irn pl emeii ting the cooling rates by 
Glover fc Jappsenl (|2007l ): in particular, we use the 
de-excitation rates calculated for the collisions with 
ionized and molecular hydrogen when available. We 
also include the cooling from neutral C, Si (as in 



Glover & J appsenl (l2007n'). Fe, and ionized O (as in 
Hollenbach fc McKed (ll98i V 

Total cooling. The cooling rate by all the above pro- 
cesses is additive and can be described by 



Atot = Ax + Ah2 + Ahd + Ago + 

-H(Ac+ -I- Ao + Asi+ + Afc+) 



(57) 



where Ax is cither Ac92 or Asd93 as appropriate, and all 
the As are functions of the temperature. Both Ah2 and 
Ahd are the cooling functions of the two molecular species, 
and the terms in brackets are the cooling of the metals. 

The overall rate of temperature change due to the heat- 
ing and cooling is given by the following equation 



dT 
dt 



7-1 



■(r-A) K/s, 



(58) 



where F is the heating source (if present), A the cooling 
source, Ui the number density (the sum is over all the el- 
ements), fcs the Boltzmann constan t, and 7 the adiabatic 
index defined in Merlin et al.l ( 201l[) as 



7 



5 IH + 5 Xb_c + 5 Xc- + 7 XHa 

3 a;H + 3 .tho + 3 x^- + 5 xn^ 



(59) 



where x is the fraction of the element indicated by the 
subscript. For an ideal gas of pure hydrogen this value 
is 5/3. If the gas is made only of molecular hydrogen we 
have 7 — 7/ 5. In eqn. ([5^ we use a linear fit of the data 
proposed bv lBolev et al.l ( 2007 ) for the adiabatic index of 
H2 under log(T) < 2.6 (see also erratum), considering a 
3:1 ortho/para ratio mix. 

The generic heating term F can be used to introduce 
heating phenomena like SNae feedback, cosmic rays, and 
others. In the current version of ROBO, F is used as a free 



parameter that is kept constant during a single run. In any 
case, the F term can be specified by the user according to 
his scientific aims. 



3. Code description 

The code has been developed with IDlH. Its user friendly 
features help the development of applications that are oth- 
erwise difficult to build in FORTRAN. The code is divided 
in self-explanatory procedures (routines) that are grouped 
in four classes (gas chemistry, dust, cooling, and general 
code behavior). The main routine calls all other subrou- 
tines that are needed to calculate the gas evolution. The 
first group of routines calculates the reaction rates and 
updates the density of each species. 

Mass conservation. The total mass per unit vol- 
ume of the ISM at time t is M{t) — X^nj(i)mi, where 
riiit) is the current number density of the i-th species, 
rrii its molecular or atomic mass, and the sum is over 
all the species. While the number densities can vary with 
time, the total mass must be conserved. In the numeri- 
cal integration of the system of differential eqns. ([T]), the 
conservation of the total mass is not always guaranteed, 
because it depends on the physical processes, the choice 
of initial parameters and the time step. Particular care 
is paid in our test computations to securing and check- 
ing run-by-run the conservation of the total mass, i.e. 
\M{to) - M{t)\/M{to) < 0(10-1°) where M{to) is the 
initial total mass of the system per unit volume. 

Time steps. Even with short time steps, it may hap- 
pen that some species reach negative values in a single 
time step. This is because the differential of a species could 
be negative, and its absolute value higher than the relevant 
number density. This problem could be solved by forcing 
the species to be positive. This would cause difficulties 
with the mass conservation, and then in the subsequent 
time step the solver would again change the sign of the 
species abundance, and finally, the cooling would produce 
some nonphysical effects. In brief, a negative value of the 
abundance of the species would artificially turn the cool- 
ing into heating, which is clearly impossible. The obvious 
way out is to suitably choose the time step. 

Numerical Solvers. In our system of differential 
equations (|T]) there are 76 reactions to deal with (64 re- 
actions plus 12 photochemical processes), so a fast and 



^ IDL is a product of ITT Visual Information Solutions, 
,http://ittvis.com/ 
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accurate solver is needed. We use the routine LSODE of 
IDL. LSODE uses adaptive numerical methods to advance 
the solution of a s ystem of ordinary differential equ ations 
by one time step (|Hindmarshlll983t IPetzoldl Il983h . opti- 
mizing the process for user-defined time steps. In our case 
LSODE absolute tolerance is 10~^°, and the relative one 
is 10-12. 

Integrators. The number densities of aU the par- 
ticles (atomic species, molecules, dusty components, and 
free electrons) are governed by the balance between cre- 
ation and destruction processes and, even more important, 
span very wide ranges of values. 

Ranges of applicability. ROBO can be safely 
used in the following intervals for temperature, density, 
and metallicity: 10 < T < 10^ K, IQ-^^ < n < 10^ and 
IQ-i^ < Z < 10. The density range is somewhat lim- 
ited towards the high value end: observational estimates 
of the density in molecular clouds can indeed reach higher 
values. Work is in progress to extend the density range. 
Furthermore, let us remind the reader that Z — 1 means 
[Fe/H] = 0, the solar case, and Z = 10 means [Fe/H] = 1; 
i.e., the ratio of the number densities (Fe/H) is equal to ten 
times solar. Finally, as not all the reaction rates cover the 
whole range of values, some extrapolations are required. 

The aforementioned ranges of applicability are chosen 
in such a way as to guarantee the numerical stability of 
the system. Indeed, the core of the model is the system of 
differential equations that describe the chemical evolution 
of the gas according to the reactions listed in Tables [5] 
through [S] The stability of the system is measured by the 
conservation of the mass of the ISM elements. Therefore, 
we carefully checked whether ROBO satisfied this con- 
ditions for values of the parameters spanning over large 
volumes of the parameter space. The results show that 
this is always the case. 

3.1. Free parameters 

ROBO contains 47 parameters governing the physics, 
the mathematics, and the numerical procedure. Here we 
briefly comment on the most important ones. 

— Ionized fraction of metals: it fixes the fraction of C, O, 
Si, and Fe at the first ionization level. 

— Metallicity: the metallicity of the gas is defined as ^ = 
X0P°/^1. It is worth recalling here that Z — 1 is the 
solar case corresponding to [Fe/H] = (see below for 
more details). 

— Number densities: the number densities of the 28 el- 
emental species and/or molecules, free electrons, and 
dusty components are all in cm"'^. 

— Fraction of carbon dust: the percentage of carbona- 
ceous grains. If it is equal to 1, all the dust is made by 
graphite grains and PAHs; if this parameter is equal to 
0, all the dust is made by silicates. Intermediate values 
are also possible. 



— Gas temperature: the temperature of the gas depends 
on the cooling and heating processes and changes dur- 
ing the gas evolution. 

— Dust temperature: the temperature of the dust is kept 
constant during each run. It controls the formation of 
the molecules on the surface of the grains and also the 
accretion of the dust grains themselves. 

— CMB temperature: the gas temperature cannot go be- 
low the temperature of the cosmic microwave back- 
ground (CMB). At present the CMB temperature is 
kept constant during each run. 

— Cosmic ray field: the field formed by the cosmic rays 
that populate the gas. It destroys or ionizes molecules 
like H2. This field is expressed in s^^, thus correspond- 
ing to the rate of ionization of the H2 by the cosmic 
rays. 

— Total integration time of a model: the total time of 
each run in s. 

— Time step: the time step used in the models. A typical 
value is lO'^ years. However, it can be changed by the 
routine LSODE. 

— Finally, there is a number of flags activating or switch- 
ing off the different sources of heating and cooling and 
the mechanisms of dust formation and destruction we 
have described in the previous sections. 

4. Models and discussion 

To validate ROBO we calculated two groups of models 
of gas evolution: the first one consists of 600 dust-free 
cases. Each case corresponds to a different combination of 
the parameters. The second group consists of 32 models 
of gas evolution, but the presence of dust grains is taken 
into account in them. 

The aim here is to investigate the evolution of the same 
elemental species in the presence or absence of the dust. 
In addition to this, the large number of cases per group 
allow us to investigate the model response at varying the 
key parameters. Indeed, owing to the large number of pa- 
rameters at our disposal, it is not possible to explore the 
whole parameter space. 

Finally, it is worth noticing that the choice of the pa- 
rameters for the two groups of models is somehow guided 
by ROBO mainly being designed to become a sort of an- 
cillary tool for NB-TSPH codes. Therefore, in view of this, 
the two groups of models use the same parameter space 
adopted by EvoL. 

4.1. Dust-free models 

Each model is a simulation of a unitary volume of gas in 
the absence of dust. The results consist of 600 (see below) 
evolutionary models of the gas temperature and number 
densities of the 28 species under consideration plus the free 
electrons. As already emphasized, in each model special 
care is paid to secure the mass conservation and to avoid 
unphysical negative number densities. The models of this 
set are calculated neglecting the presence of any type of 
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dust, and they are meant to check the internal consistency 
of ROBO. 



4.1.1. Set up of the parameters 

Each model gives the thermal and density history of the 
gas for an assigned set of parameters. The number of pa- 
rameters and the values assigned to each of them deter- 
mine the total number of models to calculate. This is given 
by = OiLi Pi I where Nm is the total number of mod- 
els. Pi the number of different values for the i-t/i parameter, 
and N the total number of the parameters, which in our 
case amounts to 47 (see Sect. 13. ip . 

For all the models, we must specify the initial value 
of the metallicity, the number densities i i ^c- j 
and the gas temperature. We adopt four values for the 
metallicity Z = 10^°/"!: {0, 10"^ IQ-^, 1}; five values for 
the H2 number density: {lO^i^ 10^6, 10-2, iQ-i, 1} cmT^- 
three values for the H+ number density: {10^^*', 10^^, 1} 
cm~^; two values for the electron number density: 
{10~^°, 10~^} cm"'^; five values for the gas temperature: 
{10, 10^, lO'', 10^ 10^} K. More details on the metallicity 
are given below when we discuss the number densities of 
the metals. All the other parameters of the models have a 
constant value. 

In all the models, the following quantities have a fixed 
initial value. They are rifj — 1-0 for the neutral hydro- 
gen, riH- = 10~^ for the hydride, — 10~^^ for the 
molecular hydrogen and nne = 0.08 for the helium. No 
helium ions are present at the beginning, so we have 
nHc+ — 'T'Ho++ = 0- The initial values of all the deuteroids 
are set to lO^^S; these are the deuterium atom, its ions 
(D+ and D^), the molecules HD and D2, and the ion HD^. 

Metals (i.e. C, O, Si, Fe and their ions C+, 0+, Si+, 
and Fe^) are computed from the metallicity as followf|f|. 
We start from the general 



and 



[Fe/H]=logior^ -logio ^ 



(60) 



where npc and riu are the number densities of iron and hy- 
drogen, and the definition of metallicity we have adopted, 
namely Z = It must be emphasized that this no- 

tation for Z is different from the commonly used definition 
of metallicity, that is, Z = 1 — X — Y with X and Y indi- 
cating the abundances by mass of hydrogen and helium. 
In the usual meaning, Z is therefore the mass fraction 
of all the species heavier than helium. In our definition 
Z is simply related to the iron content [Fe/H] by the ex- 
pression Z = 10[^°/^1. With this notation, Z can be higher 
than one: Z = 1 corresponds to the solar iron abundance. 
Considering now a generic metal indicated by X we can 
write 

(61) 



nx 



ROBO also lets us insert the value for single metallic 
species without using the total metallicity. 



nx 



Znui—] 



(62) 



Therefore if we know the ratio ('t.x/"-h)0j the total metal- 
licity Z , and the number density uh of hydrogen in the 
model, we can derive the number density of the generic 
metal X. By doing this we assume that for any metallicity 
the relative abundances of the elements follow the solar 
partition. In other words we do not consider here the pos- 
sibility that the gas may have a different distribution of 
the heavy elements from that in the Sun with respect to 
the hydrogen (a-enhancement problem). 

We use the same method to derive the number density 
of metal ions given by 



nx+ 



nx 

K Z nn — 



(63) 



where k represents the ionized fraction nx+ /nx of the 
metal X. In the present models we set k = 0, so that all the 
metals start as neutral species. If the initial temperature is 
high enough, the hot gas can ionize the metals very early 
on in the course of the ISM evolution. 

The temperature of the gas is a free parameter, but 
cannot be lower than the temperature of the CMB. For 
the present mod els we set Tqmb = 2.73 K, th e present day 
measured value ( Boeeess et al. 1992 : FixsenI 12009). 

Each model of the ISM is followed during a total 
time of 3.15 x 10^*^ s, which approximately corresponds 
to 10^ yrs. This choice stems from the following consider- 
ations. Each model of the ISM is meant to represent the 
thermal-chemical history of a unit volume of ISM whose 
initial conditions have been established at a given arbi- 
trary time and whose thermal-chemical evolution is fol- 
lowed over a time scale long enough for secular effects to 
develop but short enough to closely correspond to a sort of 
instantaneous picture of large-scale evolution of the whole 
system hosting the ISM unit volume. The initial physical 
conditions are fixed by a given set of parameters each of 
which can vary over wide ranges. One has to solve the 
network of equations for a time scale that is long enough 
to reveal the variations due to important phenomena such 
as star formation, cooling, and heating, but not too long 
to let the system depart from the instantaneous situation 
one is looking at. The value of 10^ yr resulted in a good 
compromise. 

The initial value of the time step is 3.15 x 10^° s « 
10'^ yrs. This time step determines the minimum number 
of steps required to cover the time spanned by a model. It 
means that each simulation needs at least 10"* iterations 
to be completed. The LSODE integrator may introduce 
shorter time steps depending on the complexity of the 
problem, so 10'^ is the minimum number of required steps. 
This value for the time step seems to keep the system 
stable during the numerical integration. 

While integrating, all the sources of cooling are kept 
active: r netal cooling, II 2 cooling ac cording to the p r escrip - 
tions by ' Galli fc Fallal (Il998h and IClover fc Abel (l2008h . 
and finally, the cooling from deuterated hydrogen. 
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In general, the initial values of the number densities 
fall into three groups: (i) the elemental species with con- 
stant initial values, the same for all the models (namely 
H~, HJ, He"*" and He^^); (ii) the elemental species whose 
initial values are derived from other parameters (namely 
He, all the deuteroids and the metals, such as C and O, 
which depend on the choice the total metallicity Z), and 
finally, (iii) the elemental species with free initial condi- 
tions (namely H, H2, H+, and e~). 
Hydrogen group: H, H+, H 
values of the species H^ ana ti^ are nn 



H2, and Hj . The initial 



10-9cm-3 



10- 



according to iPrieto et al. I (|2008l) . while 



the three other hydrogen species have free initial values. 
Deuterium group: D, D+, D", D2, HD, and HD+. 
The number densities of the deuteroids are calculated 
from their hydrogenoid counterparts. For the single atom 
species we have = /d "-h, "-0+ = /d JT-h+j and n^^ — 
/d"-h-i where /d = nD/".H- For the molecules, we can 
consider the ratio /d as the probability of finding an 
atom of deuterium in a population of hydrogen-deuterium 
atoms. This assumption allows us to calculate the HD, D2, 
and HD+ number densities as a joint probability. For HD 
and HD^ we have nno — /d "■H2 and nHD+ = /d '^1^+ , but 
for D2 is = /d "-112 as the probability of finding two 
deuterium atoms is /q. This is valid as long as /d <C 1. 
Helium group: He, He"*", He"*"^. The ratio nHo/nn is set to 
0.08, thus allowing the initial value of uhc to vary accord- 
ing to the initial value for njj. The initial number densities 
of the species He^, He+"'" are both set equal to zero in all 
the models. 



Metals group: C, C+, O, 0+, Si, 
Fe number density of the ISM is 



Si+, and Fe, Fe+. The 



npc = nn ■ dex [Fe/H] + log ( — 



(64) 



where (npc / n-a) q is the iron-hydrogen ratio for the Sun. 
To retrieve the number density of a given metal X we use 
nx — n-pc • fx, where fx is the metal- iron number density 
ratio in the Sun. 

The list of the species whose initial number densities 
are kept constant in all the models of the ISM is given 
in Table [TTJ The values listed here are either fixed to a 
constant value or based upon the number density of one 
of the free hydrogen species nn, riua and via the /d 
factor. Since H~ and HJ are constant, then also D~ and 
HD+ are fixed. Values are indicated as a{b) —ax 10^. 

The chemical composition of the ISM is typically pri- 
mordial with a mass abundances of hydrogen X = 0.76, 
and helium Y = 0.24 and all metals Z ~ . The helium- 
to-hydrogen number density ratio corresponding to this 
primordial by mass abundances is riHc/"-H — 0.08. The 
adopted cosmological ratio for the deuterium is /d — 
't-d/^-h — 10~^. With the aid of these numbers and the 
above prescriptions we get the number density ratios listed 
in Table 14.1.11 and the initial values of the number density 
in turn. 



Table 11. Initial values for the number densities of the 
hydrogen and helium elemental species and the deuteroids. 
See the text for details. 
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4.1.2. Results for dust-free models 

In this section we discuss the results we have obtained for 
the dust-free models. Together with this we also consider 
a side group of models calculated with some specific as- 
sumptions (e.g. without metal cooling and others to be 
explained in the course of the presentation) to underline 
some effects and to better understand the physical pro- 
cesses taking place in the ISM. 

Special attention is paid to check whether there are 
some unexpected effects or drawbacks in the models and in 
ROBO. We also describe how different physical quantities 
affect the overall behaviour of the gas. This is achieved by 
varying a single parameter at a time and keeping constant 
all the others. 

Temperature and H2 density. In Fig. 2] we compare the 
temperature evolution for different values of the initial 
molecular hydrogen densities in two cases with very dif- 
ferent metallicity, i.e. Z = lO"""^*^ and Z = 1. The important 
role played by the metal cooling is soon visible compar- 
ing the two panels. In the bottom panel (high-metallicity 
Z = 1), cooling is so strong that all the models simul- 
taneously reach in a short time interval nearly the same 
final temperature of about 10 K. In the top panel (low- 
nietallicity case) the cooling time scale is much longer than 
in the bottom panel because the contribution to cooling 
by the metals is negligible. Comparing each curve with 
the others in both panels, the effect of the molecular hy- 
drogen is evident: higher H2 densities imply steeper cool- 
ing curves and a shift toward lower temperatures. In the 
metal-free models, the temperature systematically shifts 
to lower values at increasing the density of the molecular 
hydrogen (top panel of Fig. H]). The effect is similar but 
enhanced in the case of metal-rich models (bottom panel 
of Fig. El). 

Metallicity. Figure [5] shows the evolution of the tem- 
perature (top panel) and H2 number density (bottom 
panel) of models with different metallicity Z. Each curve is 
characterized by a different initial metallicity. The sources 
of UV radiation are at work in these models. As in the pre- 
vious figure, the cooling strongly depends on the metallic- 
ity. The marked knee in the temperature-age relationship 
is caused by the cooling via the metallicity. Looking at 
the H2-age relationship, we note that the formation of H2 



T. Grassi et al.; ROBO 



21 







1 1 






\_ ^\ 














\ 1 \ - 






\ \ \ 

\ ' 




: 10° 




\ \ 




: - - 10"' 


cm"^ 


\ I 




: 10^ 


cm'^ 






: 10° 


cm'^ 






■- 10-'°Tcm-' 


1 ! 






1 1 


1 I 



log(time/yrs) 



_ 3 





1 


'1 1 ' 






























\ >\ 








\\\ 








\ 










10° 


cm'l : 






- 10"' 


cm-° : 






10-= 


cm-°] : 






10 ° 


cm-°J : 








[cmi - 




I 


■ 1 I ■ 





log(time/yrs) 

Fig. 4. Temperature evolution for different values of the 
initial H2 densities. The lines represent the time evolution 
in the cases H2 = 10"!° (solid line), H2 = 10"^ (dotted 
line), H2 = 10^2 (dashed line), H2 = 10"^ (dashed-dotted 
line), and finally, H2 = 1 (three dots-dashed line). All 
the densities are in cm~'^. Top panel: the models with 
Z ~ 10^^°. Bottom panel: the same as in the top panel 
but for Z ~ 1. See the text for more details. 



is favored when more free electrons are present (an effect 
that becomes clear when considering Figs. [51 [51 and [7]). 

UV radiation. The UV radiation plays a key role in the 
H2 density evolution. To cast light on the issue, we cal- 
culated models with no UV sources. They are displayed 
in Fig. [51 The top panel shows the temperature - age rela- 
tionship, whereas the bottom panel shows the H2 number 
density -age relationship. Comparing the bottom panel of 
Fig. [5l to that of Fig. [6l it is clear the effect of the UV 
radiation. When the UV radiation is present (Fig. [5]), H2 
is destroyed very early on, whereas when the UV radia- 
tion is absent (Fig. [5]) the number density of H2 remains 
nearly constant during the whole evolution. We also note 
that the effect of the cooling is stronger in the first case. 
This happens because the UV flux increases the density of 
free electrons, and since they act as colliders, the eventual 
higher number of collisions determines a marked improve- 
ment in the cooling process. This effect of the UV flux on 
the number density of free electrons is displayed in Fig. 
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Fig. 5. Top panel: Temperature evolution for different 
mctallicities Z and in presence of the UV radiation fleld. 
Bottom panel: evolution of the number density of molecu- 
lar hydrogen for the same set of parameters. We represent 
the cases Z = 1 (dot-dashed line), Z = lO"'^ (dashed line), 
Z = 10-2 (dotted line), and flnally Z = 10"^° (solid line). 
See the text for more details. 



[3 where we show the time evolution of this quantity at 
varying the metallicity and switching on (top panel) and 
off (bottom panel) the UV radiation. 

Ionized versus neutral metals. The difference in the 
coohng between Fig. [5l (top panel) and Fig. [51 (top panel) 
could be attributed to the different efficiencies of the cool- 
ing by ionized and neutral metals. Indeed, it is worth 
noticing that, in order to make the effects of the met- 
als evident, the initial number density of H2 in Figs. [5l 
and [His set to 10"^. The UV radiation affects the ioniza- 
tion state of atoms and molecules, thus varying the par- 
tition between ionized and neutral metals. To clarify this 
issue, we flrst compare models in the presence of UV flux 
in Fig. [51 but the metal cooling by ionized species (top 
panel) or neutral species (bottom panel) are alternatively 
switched off. The initial number density of H2 is always 
set to 10~^. For each metallicity, the values shown in Fig. 
[SI (top panel), where both neutral and ionized metals are 
included, are roughly the sum of the values presented in 
Fig. m (both panels). In the top panel of Fig. [H the cooling 
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Fig. 6. Top panel: Temperature evolution for different 
metallicities Z but in the absence of the UV radiation 
field. Bottom panel: evolution of the molecular hydrogen 
for the same set of parameters. We represent the follow- 
ing cases: Z = 1 (dot-dashed line), Z = 10~^ (dashed line), 
Z = 10-2 (dotted line), and finally, Z = 10" i° (sohd hue). 
See the text for more details. 
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Fig. 7. Temporal evolution of the electrons number den- 
sity for different metallicities Z and in presence of the UV 
radiation field (top panel, related to Fig.[S|) or without the 
UV field (bottom panel, same model of Fig. [6]). The metal- 
licities are Z = 1 (dot-dashed line), Z = 10^^ (dashed 
line), Z = 10^2 (dotted line), and finally Z = lO^^" (sohd 
line). See the text for more details. 



by ionized metals is switched off so that we would expect 
a behavior similar to that of Fig. [5] (top panel), where the 
UV radiation is switched off and the neutral metals are 
favored. It must be pointed out, however, that when the 
UV is switched off (Fig. [5]) most metals are neutral. Even 
if the ionized metals are switched off (Fig. [5]), the fraction 
of neutral metals will be much lower. Indeed, the shape of 
the curves displayed in Fig. |5] (top panel) remains similar 
to those shown in Fig. [5] (top panel), and the long cooling 
time scales of Fig. |6] cannot be reproduced. 

For the same effect of the UV flux on metals, in the 
top panel of Fig. [S] (where the UV flux is on and the 
ionized metals are off) cooling does not depend on the 
metallicity. The UV radiation easily ionizes metals like sil- 
icon, carbon, and iron: UV radiation below 13.6 eV from 
various astrophysical sources generates a UV background 
that ionizes at oms with a first ionization potential lower 
than 13.6 eV (|Maio et al.l 120071 ). Consequently, with the 
UV flux switched on, the total amount of neutral met- 
als becomes negligible. If the cooling by ionized metals 



is switched off and the contribution to the total cooling 
by neutral metals is negligible, it follows that the met- 
als almost do not contribute to cooling. When the ionized 
metals are put back and the negligible amount of neutral 
metals is dropped (Fig. |8]- bottom panel), we recover the 
situation of Fig. [S] as expected. 

Now the UV flux is switched off as in Fig. [5] and the 
metal cooling by ionized species (top panel) or neutral 
species (bottom panel) is alternatively removed. As we 
can see in Fig. |9l with no cooling by ionized metals and 
only with the contribution by neutral metals, we are able 
to reproduce the behavior observed in Fig. [5] (top panel). 
The case without neutral metals and only with cooling by 
ionized metals is not shown because it would lead to an 
almost constant temperature in all the cases, consistent 
with the picture we have just described. In brief, in this 
case, both H2 (the input number density of H2 is set to a 
lower value, thus contributing less to the cooling) and the 
ionized metals (without UV flux, as expected, most of the 
metals are neutral) little contribute to the cooling process. 
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Fig. 8. Temperature evolution for different metallicities 
Z and with the coohng contribution by ionized or neutral 
metals alternatively switched off. The UV flux is switched 
on. Top panel: no cooling by ionized metals. Bottom panel: 
no cooling by neutral metals. In both panels the meaning 
of the lines is as follows: the dot-dashed line is for Z = 1, 
the dashed line for Z = 10^^, the dotted line for Z = 10^^, 
and the solid line for Z = 10^^". See the text for more 
details. 



We can finally conclude that the partition between ionized 
and neutral metals fully explains the difference between 
Figs. [5]and|ni 

Metallicity and H2 vs. Cooling. In Fig.[TU]are plot- 
ted four different models with different combinations of 
the metallicity (very high or very low) and cooling by 
metals (included or switched off Furthermore, to ex- 
amine the effect of the H2 cooling, we set the input num- 
ber density of II2 equal to 10^"'^ cm^'^ (see Sect. l4.lTT] for 
more details on the selected input values for H2). This 
plot highlights how the cooling by metals and H2 affects 
the temperature and thus the physical status in the gas. 
The steep decrease of the temperature from the begin- 
ning of the simulations and common to all four models 
is due to the strong H2 cooling. As expected, the model 
with high-metallicity and cooling by metals (top panel of 

* In these models we switch off the C92 cooling to better 
highlight the sole effects due to different metallicities. 
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Fig. 9. Temperature evolution for different metallicities 
Z and with the cooling contribution by ionized metals 
switched off. Also, the UV flux is switched off. The mean- 
ing of the lines is as follows: the dot-dashed line is for 
Z = 1, the dashed line for Z = 10~^, the dotted line for 
Z = 10-2, and the solid line for Z = 10-^°. See the text 
for more details. 



Fig. [TUl dashed line) is the one experiencing the strongest 
cooling. The two models with very low-metallicity have a 
similar behavior independent of the presence or absence of 
the cooling by metals. The last case, with high-metallicity 
and no cooling by metals, is the one with the lowest tem- 
perature decrease. The reason for it is explained in Fig. 
[TTl where the evolution of H2 for the four cases under 
consideration is shown. 

The model with high-metallicity and no cooling by 
metals has the lowest H2 density compared to the other 
three and consequently has low cooling. From Fig. 1111 we 
also see that the model with high-metallicity and cooling 
by metals (dashed line) has the highest II2 density. This 
can explain the results in the top panel of Fig. [TUl in the 
sense that the cooling process here could be mainly due 
to the contribution coming from the molecular hydrogen. 
The bottom panel of Fig. [TU] shows the same models with 
no cooling by H2. In this case only the model with cooling 
by metals included and high-metallicity undergo a signif- 
icant cooling. Indeed, if the effect of H2 is neglected, we 
need metals in significant amounts to have strong cooling. 

4.2. Models with dust: parameter set up 

In this section we consider models in which the effect of 
the dust grains is taken into account. We adopt here the 
same values for the parameters as in the dust-free models. 

The minimum initial metallicity is Z ^ 10"^, and the 
initial number density of the molecular hydrogen is — 
10^® cm^^. The number densities of electrons and H+ 
are free parameters; they are nH+ = {10^^", 10^^} and 
n^- — {10-^°, lO"^}, both in units of cm^"^ as usual. The 
initial temperature is set to T = 10^ K. 
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Fig. 10. Temperature evolution with different metallici- 
ties Z and cooling options by metals. The meaning of the 
lines is as follows: high Z and no cooling by metals (solid 
line), low Z and no cooling by metals (dotted line), high Z 
and cooling by metals included (dashed line), low Z ^ and 
cooling by metals included (dashed-dotted). The dashed- 
dotted and dotted lines overlap. Top panel: the H2 cooling 
is enabled. Bot tom panel: the H2 cooling is switched off. 
The cooling by [Ce3 (I1992I) is switched off in both panels. 



The only difference with respect to the previous mod- 
els is the presence of dust. We adopt four values for 
the number density of dust grains, namely njust = 
{0, 10"'^, 10~^, 10^^} cm^-^. In these models the compo- 
sition of the dust mixture is 50% carbonaceous grains and 
50% sihcates. 

We also calculate models with or without dust sput- 
tering by shocks, in order to describe the behavior of a 
turbulent gas particle with and without the grain deple- 
tion due to the shocks. 

All the other parameters remain the same as in the 
dust-free models, such as the number densities of different 
elements and the cooling processes. In all these models 
thermal sputtering and dust formation are active, so the 
dust properties are let change during the evolution of the 
interstellar medium. 
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Fig. 11. Evolution of the number density of molecular hy- 
drogen II2 for different initial metallicity Z and options for 
the cooling by metals. The meaning of the lines is as fol- 
lows: high Z and no cooling by metals (solid line), low Z 
and no cooling by metals (dotted line), high Z and cool- 
ing by metals enabled (dashed line), low Z and cooling 
by metals switched off (dashed-dotted line). The dashed - 
dotted and dotted lines overlap. The cooling bv lCen (1992) 
is switched off in both panels. 



4.2.1. Results for models with dust 

The series of plots going from Fig.lT^through Fig.lTilshow 
the temporal evolution of three important quantities of the 
models, namely the gas temperature, the number densi- 
ties of dust, and molecular hydrogen for different amounts 
of initial dust and different initial metallicities. Clearly 
there is a tight relationship between the initial amount 
of dust, the temporal behavior of temperature, and the 
number density of H2. First, at a given low-metallicity 
and by increasing the dust fraction, the temperature de- 
creases earlier and faster (top panel of Fig. [T^ . whereas 
if the metallicity is high there is no remarkable effect of 
the increased dust content (bottom panel of Fig. fT^ . Also, 
for the low-metallicity, the trend is anticipated as the dust 
content increases. 

Looking at the temporal evolution of the number den- 
sity of H2 shown in the panels of Fig. [T3] we note that, in 
coincidence with the temperature fall off and subsequent 
gentle decrease, the H2 number density first decreases and 
then increases, forming a local minimum. This happens 
because the dust drives the formation of H2 : models with 
more dust suffer more cooling and hence form more H2. 
Only if the metallicity is very high, does the cooling effect 
of dust grains lose importance as shown by the bottom 
panel of Fig. [T^l To conclude, the curves displayed in the 
top and bottom panels of Fig. [13] nearly have the same 
shape, but a different vertical offset. The position and 
depth of the minima depends on the temperature varia- 
tion, and the vertical offset clearly depends on the amount 
of dust present in the gas. This means that amount of dust 
and H2 are closely related. 
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Fig. 12. Temperature evolution for different amounts of 
initial dust. Top panel: low-metallicity cases. Bottom 
panel: high metallicity cases. In both panels the meaning 
of the lines is as follows: the solid line is for rtdust — 10^'^, 
dotted line is for Udust = 10~^, and finally, the dashed line 
is for ndust = 10~^. See the text for more details. 



Finally, we calculate and present six models (three for 
each initial amount of the dust density) that include now 
the dust destruction by shock sputtering and vaporiza- 
tion. Their temporal evolution is limited to the first 10^ 
years only. The results are shown in Fig. [15] (the gas tem- 
perature) and Fig. (the dust grains number density as 
a function of the size of dust grains) . Looking at Fig. [151 
the temporal evolution of the gas temperature is the same 
as in the previous case with the sole thermal sputtering at 
work. In contrast, the dust number density undergoes big 
changes as shown in Fig. [THl because the shock sputtering 
is very efficient. Figure [THj shows how the dust distribution 
changes with the time for the rtdust = 10~^ cm~^ case. 
After plotting the data of Fig. [THl the dust grains have 
been grouped in bins according to their sizqfj. The lines 
represent different distributions at different ages, namely 



® It worth noticing that even if the initial distribution of dust 
density per size bin is a power law, assuming that the dimension 
of the bin size is chosen in such a way that the density per bin 
is constant, the initial distribution shown in Fig. [16] appears to 
be flat. 
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Fig. 13. Evolution of molecular hydrogen for different 
amounts of initial dust. Top panel: low-metallicity cases. 
Bottom panel: high metallicity cases. In both panels the 
meaning of the lines is the same as in Fig. [121 See the text 
for more details. 



10, 10^, 10^, and lO'* years from top to bottom. As ex- 
pected, the shock first destroys the large size grains and 
then the small ones. 

5. Including ROBO results in an NB-TSPH code 

In this section we briefly describe how ROBO is planned 
to be used as an ancillary tool for the NB-TSPH code 
EvoL to calculate the thermal and chemical properties of 
the gas particles. In the following we present some prelim- 
inary results. The method and it s results will be w idely 
discussed in a forthcoming paper ( Grassi et al. 201 1[ ). 

We already mentioned that two different techniques 
can be used. The first one is a real-time method, in 
which the elemental abundances are calculated solving the 
Cauchy problem of the chemical network. It means that 
for every time step and for every particle, EvoL needs to 
perform this computation using ROBO as a routine ded- 
icated to this purpose. The advantage is that elemental 
abundances, thermal properties of the gas particles, and 
cooling can be calculated with high precision. However, 
this would require large computational resources, and un- 
fortunately, the more elements it tracks, the more CPU 
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Fig. 14. Dust density evolution for different initial 
amounts of dust. Top panel: low-metallicity cases. Bottom 
panel: high metallicity cases. The meaning of the lines is 
the same as in Fig. [121 See the text for more details. 
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Fig. 15. Temperature evolution for different amounts 
of dust. Shock disruption is enabled. Top panel: lo'w- 
metallicity cases. Bottom panel: high-metallicity cases. 
The meaning of the lines is the same as in Fig. [12] See 
the text for more details. 



time is needed. Furthermore, the net"work complexity and 
additional physical processes (e.g. detailed cooling) would 
increase the computational cost of it. 

To avoid these problems an NB-TSPH code can use the 
s- called "grid" method. It consists of calculating a large 
number of models exploring different scenarios produced 
with different sets of the input parameters. The results 
are then tabulated. When a particle in the cosmological 
or galaxy simulation needs to know its evolution after a 
time step, the particle looks for the closest parameters set 
in the grid, and interpolate the evolved parameters. This 
a lightweight method, but the main drawback is that the 
number of free parameters cannot be too large, mainly be- 
cause interpolating multidimensional grid can be difficult 
and inaccurate (depending on the grid coarseness). 

Our method is similar to the latter one. First of all, we 
select a set of parameters. Second, we run ROBO for each 
parameter combination over the total integration time we 
have chosen (3.15 x 10^"' s). The total number of simula- 
tions depends on the number of parameters and on the 
multiplicity of each parameter (see Sect. I4.1TT|) . This part 
can easily be achieved using for example a cluster of com- 
puters. 



Now we have a huge number of models characterized 
by a large set of parameters. With the standard grid 
method, interpolating the requested data would be com- 
plicated. To solve this problem we use a powerful tool: 
the ANNs. There are different kinds of ANN architectures 
depending on the task it needs to perform. We c hoose the 
so ca lled back-propagation algorithm (Rumelhart et al. 
19861) . It is one of the most used and versatile methods 
for solving a wide range of problems. This mathemati- 
cal method consists of an algorithm that behaves like the 
gas model itself. Our ANN learns how the gas evolves us- 
ing the models produced with ROBO. We first perform 
a training stage, in which a subset of randomly chosen 
models is shown to the ANN. After some iterations the 
ANN converges to a stable state. At this point we look at 
the complementary of the subset of models: its aim is to 
calculate the ANN error on unknown data. If the error is 
small enough, we can assert that the ANN behaves like 
the gas model. When this happens, in EvoL ROBO can 
be safely replaced by the ANN. 

The great advantage is that the ANN can give the re- 
quested solution with ease. It is very light, since it consists 
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Fig. 16. Dust density evolution for different dimension 
bins. Shock disruption is enabled. Lines represent the dis- 
tribution of dust among the bins for different ages (from 
top to bottom, 10, 10^, 10'^, and 10'* years). See the text 
for more details. 
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Fig. 17. Comparison between the theoretical data calcu- 
lated by ROBO for the free electron number density and 
fed to the ANNs (squares) and the same data retrieved 
by the ANNs (crosses). Each input data corresponds to a 
model calculated by ROBO for a set of input parameters, 
hence to a network data point of the hyperspace of param- 
eters handled by the ANNs. The input and predicted data 
refer to the free electron number density after 10^ years 
of evolution. The data are normalized to the maximum 
value of the distribution. Similar plots are easily obtained 
for all the other physical quantities relative to the ISM. 



of three small matrices of so-called synaptic weights. Each 
particle needs only to perform a matrix product, which is 
(in this scenario) a fast algebraic operation. In this way 
retrieving chemical evolution for each particle is fast. To 
illustrate the ability of the ANN to provide good repre- 
sentations of the physical state of the ISM in Fig. [T71 we 
show the comparison between the theoretical data calcu- 
lated by ROBO for the free electron number density and 
fed to the ANN (squares) and the same data retrieved by 



the ANN (crosses). The agreement is remarkably good. 
The same technique can be applied to any other variable 
of interest here. 

6. Conclusions 

We have presented a model of the ISM that provides a de- 
tailed description of the gas chemistry and evolution, the 
formation and destruction of dust grains of different types, 
and finally, a thorough description of the cooling process 
over a wide range of physical parameters and initial con- 
ditions. The way the model is conceived corresponds to 
an instantaneous picture of the physical state of an el- 
ementary volume of the ISM characterized by a set of 
physical parameters assumed here as the initial conditions 
of a given volume element. Under the action of the ISM 
models, the initial physical state evolves on a secular time 
scale. This provides us a sort of vector field telling how a 
given physical state will evolve (how much and in which 
direction in the multidimensional space of the physical 
conditions). The integral of the elementary volumes over 
the underlying evolutionary path of the grand physical 
quantities like density and temperature (all of these func- 
tions of space and time) of the host system (a galaxy or 
a cosmological simulation) will give us the detailed evolu- 
tion of the ISM. This is the big advantage offered by the 
model, securing it a wide range of applicability. 

We have presented here the results for dust-free and 
dust-rich ISM at varying the key parameters. The first 
group of models for a dust-free medium is meant to under- 
stand how the ISM behaves in the absence of dust grains. 
These models highlight the importance of the different 
kinds of cooling that are dominant in different kinds of 
environments. In particular, we call attention to the role 
of metals and free electrons in driving the physical behav- 
ior of the ISM via their effect on the gas cooling during 
its evolution. The dust-rich ISM allows us to understand 
how the ISM responds to the presence of the dust. In par- 
ticular, we analyzed the temperature variations caused by 
the presence of dust in different amounts. We have also ex- 
plored how the creation and destruction of the dust grains 
(the latter induced by shock and thermal sputtering) af- 
fects the evolution of the ISM. 

The ISM model and companion code were created as 
auxiliary tools for NB-TSPH simulations in the context 
of galaxy cosmological simulations of the Universe and 
models of galaxy formation, structure, and evolution. Our 
specific aim is to give a more accurate description of the 
gas component in EvoL or in similar codes in litera- 
ture. Finally, ROBO is also designed to run in small and 
middle-size computers. Thanks to ROBO, detailed gas 
physics can be inserted in NB-TSPH simulations at low 
computational costs. 

To include the resuhs of ROBO in our NB-TSPH 
code, we plan to use the ANNs that are more accurate, 
faster, and easier to implement than the standard fits on 
multidimensional g rids. A complete a ccount of this will be 
made public soon (jGrassi et al.ll201l[ ). 
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Future implementations of ROBO are planned, among 
which we mention the inclusion of the photo-ionization 
by single stellar populations of different age and chemical 
composition in the chemical network and a better deter- 
mination of the grains temperature that is tightly related 
to the local stellar radiation field. 
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